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ABSTRACT 

The  vertical  discretization  in  a  linearized  baroclinic  prediction  model  was  analyzed 
by  comparing  various  finite  element  and  finite  difierence  solutions  following  Jordan 
(1985)  and  Shapiro  (1987).  The  baroclinic  instability  experiments  of  Shapiro  (1987)  were 
augmented  to  include  the  unstaggered  vertical  scheme  from  Jordan  (1985).  Two  basic 
wind  profiles  were  used  and  the  experiments  were  run  with  different  resolution  models 
and  horizontal  wavelengths.  For  a  given  wind  profile  and  vertical  resolution,  different 
models  performed  better.  The  finite  element  models  for  the  staggered  vertical  grids  did 
not  perform  up  to  their  possibiUties  due  to  the  boundary  elements.  However,  for  the 
unstaggered  vertical  grid,  the  finite  element  model  did  better  than  the  finite  difierence 
model  in  most  cases. 
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I.     INTRODUCTION 

At  the  present  time,  most  numerical  weather  prediction  models  use  fmite  differences 
to  accomplish  the  vertical  discretization  even  though  they  use  fmite  difference,  fmite  el- 
ement, or  spectral  horizontal  discretizations.  The  Canadian  regional  and  hemispheric 
models  (Staniforth  and  Daley,  1977  and  1979)  which  use  fmite  elements  in  the  vertical 
are  an  exception.  The  successful  numerical  prediction  of  synoptic  evolutions  requires  a 
proper  representation  of  the  vertical  variation  of  the  predictive  fields.  Since  smaller  scale 
features  such  as  fronts  (Hoskins  and  Bretherton,  1972  and  Williams,  1967)  and  the  large 
scale  planetar}'  waves  (Gall,  1976)  are  forced  by  energetic  synoptic-scale  features,  it  fol- 
lows that  all  predictive  scales  of  motion  may  be  sensitive  to  the  vertical  discretization 
used  in  the  numerical  mdels. 

Most  of  the  fmite  difference  vertical  discretizations  use  a  staggered  arrangement  of 
variables.  Winninghoff  (196S),  Arakawa  and  Lamb  (1977)  and  Schoenstadt  (1980) 
demonstrated  that  staggering  of  variables  in  the  horizontal  improves  geostrophic  ad- 
justment and  the  response  to  small  scale  forcing.  Most  quasi-geostrophic  models 
(Charney  and  Phillips,  1953)  use  vertical  staggering  where  the  vertical  motion  and  the 
temperature  are  carried  between  the  levels  which  carry  horizontal  velocity  and  pressure. 
This  arrangement  will  be  referred  to  as  grid  B.  Lorenz  (1960)  introduced  a  different  grid 
for  the  balance  equations  which  was  designed  to  conserve  energy.  This  arrangement 
places  only  the  vertical  velocity  between  the  levels  which  carr\'  the  other  variables  (hor- 
izontal velocity,  pressure  and  temperature)  and  will  be  referred  to  as  grid  A.  Tokioka 
(1978)  analyzed  a  number  of  vertical  grids  with  linearized  equations  and  found  that  grid 
A  has  a  computational  mode  in  the  temperature  field.  Arakawa  (1984)  compared 
baroclinic  instability  for  grids  A  and  B  in  the  linearized  quasi-geostrophic  equations. 
He  found  a  false  short  wave  instability  for  grid  A  which  did  not  occur  with  grid  B.  This 
problem  is  related  to  the  computational  mode  in  the  temperature  field.  Another  diffi- 
culty with  grid  B  is  that  the  matrix  which  must  be  inverted  to  fmd  the  temperature  from 
the  pressure  is  singular.  This  is  especially  important  for  initialization.  Many  operational 
primitive  equation  models  use  grid  A  for  energy  conservation  (Shapiro,  1987). 

The  use  of  finite  elements  for  the  vertical  discretization  can  be  expected  to  give  a 
more  accurate  representation  of  vertical  variations.  The  finite  element  method  is  a  spe- 
cial case  of  the  Galerkin  procedure  which  represents  the  dependent  variables  with  a 


weighted  sum  of  basis  functions  that  have  a  prescribed  spatial  structure.  The  finite  ele- 
ment method  employs  basis  functions  which  are  zero  except  in  a  limited  region  where 
they  are  low-order  polynomials.  This  method  has  been  used  in  engineering  statics  (e.g., 
Zienkiewicz,  1977)  and  it  has  been  applied  to  fluid  dynamics  and  hydrology  (see  Gray 
and  Pinder,  1976).  The  finite  element  method  has  been  successfully  apphed  to 
meteorological  prediction  with  the  shallow  water  equations  by  Cullen  (1973),  Hinsman 
(1975)  and  Staniforth  and  Mitchell  (1977,  1978).  Cullen  (1973),  Neta  et  a/  (1986),  and 
Neta  and  Williams  (1986)  demonstrated  that  finite  element  formulations  with  piecewise 
linear  basis  functions  are  more  accurate  than  second  order  finite  difierences. 

Jordan  (1985)  compared  six  linear,  baroclinic,  vorticity-divergence  equation  models 
using  three  grid  schemes,  grid  A,  grid  B  and  an  unstaggered  grid.  The  three  grids  are 
depicted  in  Fig.  1.  Shapiro  (1987)  improved  the  models  for  grids  A  and  B,  and  corrected 
several  problems  in  the  heating  and  mountain  forcing  experiments.  He  also  ran  a 
baroclinic  instability  experiment  without  heating  or  mountain  forcing  using  linear  verti- 
cal shear  in  the  mean  wind  profile.  With  the  ever-increasing  density  of  observations  and 
the  advent  of  satellite-based  cloud  pictures,  transient  disturbances  of  subsynoptic  scales 
have  been  observed  in  certain  preferred  locations  such  as  polar  lows  (Satyamurty  et  al, 
1982).  Satyamurty  et  al  (1982)  found  that  unstable  modes  can  be  generated  at  smaller 
horizontal  wavelengths  using  hyperbolic  tangent  profiles.  Farrell  (1985)  suggests  that 
for  backward  tilting  initial  states  the  waves  can  experience  rapid  growth  in  the  early 
stages  and  then  settle  into  the  most  unstable  growth. 

The  purpose  of  this  study  is  to  augment  the  barochnic  instability  study  by  Shapiro 
(1987)  by  including  the  unstaggered  models  from  Jordan  (1985)  and  by  looking  at 
smaller  resolution  models  and  smaller  wavelengths.  Various  hyperbolic  tangent  vertical 
profiles  of  the  mean  wind  field  are  used  to  focus  on  more  difiicult  baroclinic  instabiUty 
situations  with  smaller  horizontal  wavelengths  and  vertical  scales,  and  with  greater  sen- 
sitivity to  large  local  shears.  This  study  also  experiments  with  initial  states  which  have 
a  backward  tilt  with  height  in  an  effort  to  investigate  the  enhanced  initial  growth  due  to 
the  tilt.   The  results  of  the  experiments  are  given  in  Chapter  III. 
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Fig.   1.      Three  vertical  grids. 


II.     MODEL  DESCRIPTIONS 

A.  MODEL  FEATURES 

Jordan  (1985)  developed  six  numerical  models  with  several  features  to  make  easy 
modifications  for  a  wide  range  of  experiments.  Shapiro  (1987)  modified  the  A  and  B 
models  to  improve  the  performance  of  the  models  near  the  boundary  and  used  numerical 
integration  to  evaluate  the  entries  of  the  matrices.  He  also  added  menus  to  make  the 
transition  between  four  experiments  simple.  The  user  is  able  to  prescribe  heating, 
mountain  topography,  velocity  perturbation,  or  baroclinic  experiments  and  the  model 
will  make  the  prescribed  changes  in  the  variables  governing  these  cases.  The  vertical 
velocity  vector  in  finite  element  model  B  has  been  corrected  to  include  an  additional  el- 
ement in  an  attempt  to  eliminate  the  oscillation  in  the  temperature  profile  from  Shapiro 
(1987).  The  two  models  for  grid  C  have  been  modified  to  run  baroclinic  instability  ex- 
periments and  to  include  the  same  menus  as  the  A  and  B  models.  Another  menu  has 
been  added  to  all  six  models  to  select  between  linear  and  nonlinear  wind  shear  profiles. 
The  models  are  written  in  modular  structure  using  FORTRAN  77,  There  is  parallel 
construction  between  models.  The  subroutines  used  in  one  model  are  very  similar  to 
those  used  in  the  other  models.  The  models  can  run  quickly  on  an  IBM-3033 
mainframe.  For  example,  a  96-hour  forecast  for  a  6-layer  fmite  element  model  uses  less 
than  six  seconds  of  computer  processing  time;  however,  the  same  forecast  for  60  layers 
can  take  up  to  220  seconds. 

B.  GOVERNING  EQUATIONS 

Each  model  approximates  the  same  set  of  governing  equations.  The  vorticity 
equation  (2.1),  the  divergence  equation  (2.2)-,  the  surface  geopotential  equation  (2.3), 
and  the  first  law  of  thermodynamics  (2.4)  are  the  prognostic  equations  for  the  forecast 
variables  vorticity,  divergence,  surface  geopotential  and  potential  temperature.  The 
surface  geopotential  equation  is  the  lower  boundary  condition  on  the  vertical  velocity. 
The  vertical  coordinate  Z  =  —  ln(p/po)  is  used,  but  the  non-boussinesq  terms  involving 
e-^  are  replaced  by  one.   The  prognostic  equations  in  the  coordinates  x,  y,  Z  and  t  are 

-^  +  (C  +  0D  +  /iv  +  i^-|^-42^-?^  =  0,  (2.1) 

dt  ex    oZ        oy    oZ 


dD    ,  ,  cu  .2  ,   ,  r\'2  ,   ^  cu    dv    ,    cw    cu    ,    dw    c\ 


ex  c\  c\    ox 


+  — — ^rr  + 


dt  ex  cy  dy    dx        dx    dZ        cy    rZ 

+  /?u  -  f^  +  V^(/)  =  0,  (2.2) 

'  =MTSand  (2.3) 


dt 
dT 


dt 


=  Q-  (2.4) 


In  these  equations, 

C  is  the  vertical  component  of  vorticitv,  C  =  — r~  , 

ox       cy 

cu       cv 
D  is  the  horizontal  divergence,  D  =  — — h  -::;—  , 

dx       oy 

4>  is  the  geopotential,  <^  =  gZ, 

</)s  is  the  surface  geopotential, 

T  is  the  potential  temperature, 

u  is  the  x-component  of  velocity, 

V  is  the  y-component  of  velocity, 

w  is  the  vertical  velocitv, 

Q  is  the  diabatic  heating  per  unit  time  per  unit  mass, 

MTS  is  the  forced  vertical  velocity  due  to  flow  over  mountain  topography, 

f  is  the  Coriolis  parameter, 

j5  is  df  dy, 

d()       d{)  d{)  c{)  d{) 

+  u  -T v  V  —: —  +  w  —r^r      and 


dt  dl  dx  dy  dZ 

V^  is  the  horizontal  Laplacian  operator. 

The  prognostic  equations  are  linearized  by  expanding  the  variables  into  their  mean 
and  perturbation  states,  as  in  Jordan  (1985).  The  resulting  linearized  forecast  equations 
are 

dV  dV 

-T-  =  -fD'-u^r--/?v',  (2.5) 

ot  ox  ^ 


dP'  ^.      -  cD'       a  ,       dw'    du        o  0                                                              f^,. 

ct  ex                   ^x    dZ        ^x 

8(b/  d<bJ            d(h^         — 

'  =  -  u  -T^  -  V'  -^  -  RTw'  +  MTS'  and                                                    (2.7) 


^t  ex  dy 

rT'  rT'  rT  rJ 

-ar=-"-ir-^iF-"iz+"2.  (2-8) 

where  R  is  the  gas  constant  for  air,  (')  denotes  perturbation  quantities  and  (  )  denotes 
mean  quantities.  The  use  of  Xbar  in  the  text  will  be  used  to  denote  mean  quantities  of 
a  variable  X  (Shapiro,  1987). 

The  diagnostic  variables,  u',  v',  w'  and  C',  are  calculated  from  the  forecast  variables 
using  the  definitions  of  divergence,  vorticity,  the  hydrostatic  equation  and  the  continuity 
equation.   The  relationships  are  given  in  equations  (2.9)  through  (2.12). 

^u' 


=  D'.  (2.9) 

ox 


(2.10) 


d\ 


e<p' 

-T^=RT'.  (2.11) 

cZ 

D'+4^  =  0.  (2.12) 

The  use  of  primes  to  denote  perturbation  quantities  will  be  discontinued.  All  quantities 
used  in  the  remainder  of  the  paper  will  be  perturbation  quantities  unless  otherwise  noted. 
The  mean  state  is  assumed  to  be  in  hydrostatic  and  geostrophic  balance.  The  term 
oTbar/^y  in  the  first  law  of  thermodynamics  can  be  evaluated  by  taking  d!dy  of  the 
hydrostatic  equation  and  substituting  for  50bar/5y  from  the  geostrophic  relation, 
<9</)bar/^y  =  —  f ubar.   Thus, 

-^  =  -  — — .  (213) 

dy  R   cZ  ^       ^ 

Geostrophic  balance  of  the  mean  state  at  the  surface  implies 


-T-  =  -^^sfc-  (2.14) 

The  expressions  (2.13)  and  (2.14)  are  substituted  into  equations  (2.8)  and  (2.7),  respec- 
tively. 

A  singlewave  spectral  representation  is  used  in  the  x-direction,  with  wave  number 
/!=  2;r;L,  where  L  is  the  wavelength  in  the  x-direction.  The  perturbation  quantities  have 
the  form 

C(x,Z,t)  =  Aj(Z,t)  cos  /ux  -I-  A2(Z,t)  sin  ^^x,  (2.15) 

D(x,Z,t)  =  Di(Z,t)  cos  ^x  +  D2(Z,t)  sin  ^x,  (2. 1 6) 

T(x,Z,t)  =  Tj(Z,t)  cos  MX  +  T2(Z,t)  sin  ^x,  (2.17) 

(^s(x,Z,t)  =  S,(Z,t)  cos  A^x  +  S2(Z,t)  sin  ^x,  (2.18) 

u(x,Z,t)  =  Ui(Z,t)  cos  ux  +  U2(Z,t)  sin  mx,  (2.19) 

v(x,Z,t)  =  V,(Z,t)  cos  ^x  +  V2(Z,t)  sin  fix,  (2.20) 

w(x,Z,t)  =  \\\{Z,i)  cos  AiX  +  W2(Z,t)  sin  fxx,  (2.21) 

0(x,Z,t)  =  H|(Z,t)  cos  fix  +  H2(Z,t)  sin  /^x,  (2.22) 

Q(x,Z,t)  =  Qi(Z,t)  cos  pix  +  Q2(Z,t)  sin  ^x  and  (2.23) 

MTS(x,Z,t)  =  MTSi(Z,t)  cos  lux  +  MTS2(Z,t)  sin  ^x.  (2.24) 

The  relations  (2.15)  through  (2.24)  are  substituted  into  equations  (2.5)  through 
(2.12).  The  prognostic  and  diagnostic  equations  are  separated  into  equations  for  the 
cosine  and  the  sine  terms.   The  resultant  prognostic  equations  are 

-^  =  -fD,  -  u/xA2  -  ^V„  (2.25) 

-^  =  -nD2  +  uM,-A^2.  (2.26) 


^D]       ,,        _  ^        ^,.  du 
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^^    =fA, -UA.D2-/?L\-/z-^W2  +  /Hi,  (2.27) 

-^  =  fA2  +  U;uD,  -  ;5U2  +  /^  -^  Wi  +  m'H2,  (2.28) 

4r^  =  -  u/^Sj  +  f  U3f,V,  -  RTW,  +  MTS,  and  (2.31) 

dS  — 

-TT-  =  umS]  +  f  Usf,V2  -  RTW2  +  MTS2.  (2.32) 

The  resultant  diagnostic  equations  for  u  and  v  are 

U,=--7r,  (2.33) 

^2  =  11-'  (2-3'4) 

A. 

Vi  =  --^   and  (2.35) 

V2  =  4^.  (2.36) 

Geopotential  values  above  the  surface  are  obtained  by  integrating  the  hydrostatic 
equation  from  the  surface  (Z  =  Zq)  to  height  Z. 


H,  =  R|    Ti(Z,t)dZ  +  Si  and  (2.37) 


H2  =  rJ  T2(Z,t)dZ  +  S2.  (2.38) 


The  vertical  velocity  is  calculated  by  integrating  the  continuity  equation  from  the  top 
of  the  atmosphere  (Z  =  Zj)  down  to  height  Z.  The  upper  boundary'  condition,  \v  =  0 
at  Z  =  Zt,  is  used.   The  diagnostic  equations  for  the  vertical  velocity  are 

W)  =  I    ^Di(Z,t)dZ  and  (2.39) 


W2  =  J  ^D2(Z,t)dZ.  (2.40) 


Equations  (2.25)  through  (2.40)  are  the  prognostic  and  diagnostic  equations  that 
govern  all  six  numerical  models.  Using  the  given  basic  state  and  the  one-wave  spectral 
perturbation  quantities,  the  governing  equations  reduce  to  functions  of  Z  and  t.  The 
models  are  effectively  one-dimensional  (Shapiro,  1987). 

To  display  the  results  of  each  model,  the  sine  and  cosine  amplitudes  of  each  variable 
are  combined  to  determine  the  amplitude  and  phase  of  a  single  cosine  wave  in  the  x- 
direction.  A  typical  variable  has  the  form: 

Y(x,Z,t)  =  A(Z,t)  cos(mx  -  ^),  (2.41) 

where  the  amplitude  is  A(Z,t)  and  the  phase  is  (5(Z,t).  The  amplitude  and  phase  are 
calculated  at  each  level  for  all  variables. 

C.     TIME  DIFFERENCING 

Tv^'o  forward  time  steps  are  taken  to  start  each  model  and  then  leapfrog  time  dif- 
ferencing is  used.  The  leapfrog  scheme  is  a  higher  order  scheme  than  forward  difTer- 
encing  and  is  the  one  used  in  many  numerical  models.  A  Robert  filter  is  used  to  reduce 
the  amplitude  of  the  computational  mode  generated  by  the  leapfrog  time  differencing. 
The  filter  is  discussed  by  Haltiner  and  Williams  (1980).  For  a  prognostic  variable  F, 
calculate  Fbar^.i,  the  average  value  of  F  at  time  step  (n-l)At,  using  equation  (2.42), 

F„_,  =  F,_,  +  y(F,  -  2F„_,  +  F^.^),  (2.42) 

where  y  is  a  weighting  function.  Using  the  unaveraged  values  at  time  step  nAt,  compute 
the  tendency  {d¥ldi)^  from  its  predictive  equation.  The  predicted  value  at  time  step 
(n+  l)At  is  then  calculated  using  equation  (2.43), 


F,^,  =  F,_,  +  2At(^)„.  (2.43) 

In  all  the  experiments,  y  =  0.05  is  used.  The  time  step  for  each  experiment  is  calculated 
in  the  model  by  requiring,  for  computational  stability, 

vAt  =  -^,  (2.44) 

where  v  =  /xc  and  c  is  the  typical  phase  speed  of  an  external  gravity  wave. 

D.  VERTICAL  GRIDS 

Each  of  the  models  uses  one  of  three  vertical  grids.  The  three  methods  for  distrib- 
uting the  variables  over  discrete  levels  are  depicted  in  Fig.  1.  The  staggered  levels  are 
represented  by  the  dashed  lines  in  Fig.  1.  Notice  that  the  heights  at  which  the  variables 
are  defined  change  between  the  three  grids.  The  notation  used  in  this  paper  to  denote 
the  staggered  and  unstaggered  levels  is  consistent  with  the  conventions  used  in  the  coded 
models.  The  height  of  the  unstaggered  levels  is  denoted  as  Z'.  The  height  of  the  stag- 
gered levels  is  denoted  as  Z.  In  the  models,  both  Z/  and  Z,  are  defined  to  be  the  surface 
of  the  earth.  It  is  assumed  that  the  staggered  level  Z,  is  exactly  in  the  middle  of  the  layer 
between  Z,_,'  and  Z,'.  This  distinction  is  important  because  the  models  can  have  layers 
with  unequal  depth.  Thus,  the  height  of  the  staggered  levels  is  defmed  relative  to  the 
height  of  the  unstaggered  levels. 

A  finite  diflference  model  is  written  for  each  of  the  grid  structures.  The  models  are 
denoted  as  FDM-A,  FDM-B  and  FDM-C.  Similarly,  finite  element  models  using  the 
three  grids  are  indicated  by  FEM-A,  FEM-B  and  FEM-C. 

E.  FINITE  DIFFERENCE  MODELS 

The  only  differences  in  the  equations  between  the  three  FDM  models  are  the  ap- 
proximations of  terms  involving  dubar/dZ  and  5Tbar/5Z  in  the  prognostic  equations  and 
the  approximations  of  the  integral  in  the  diagnostic  geopotential  equation.  Centered 
difference  approximations  are  used,  except  at  the  boundaries  where  one-sided  differences 
are  used.  The  finite  difference  approximations  used  in  the  prognostic  equations  are 
listed  in  Appendix  A. 
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F.     FINITE  ELEMENT  MODELS 
L     FEM-A 

The  FEM-A  model  defines  vertical  velocity  (w)  at  the  unstaggered  levels  in 
terms  of  the  basis  functions  lAXZ).  The  other  variables  are  defined  at  the  staggered  levels 
in  terms  of  the  basis  functions  4),{Z).   The  expansion  for  a  typical  term  is 


n+l 

Ai(Z.t)  =  N  Ai(t).Aj(Z)       or 


j=i 


n+l 


1 


J=l 


A\mm. 


(2.45) 


The  basis  functions  for  this  model  are  depicted  in  Fig.  2.  The  basis  functions  \l/(Z)  are 
defined  for  the  unstaggered  levels  (solid  lines  at  height  Z')  and  the  basis  functions  4>{Z) 
are  defined  for  the  staggered  levels  (dashed  lines  at  height  Z). 


Fig.  2.      Basis  functions  for  grids  A  and  B. 

The  finite  element  approximations  for  the  vorticity,  divergence  and 
thermodynamic  equations  are  derived  by  substituting  the  expansion  for  each  dependent 
variable  into  equations  (2.25)  through  (2.30).    Each  equation  is  multiplied  by  </',(Z)  and 
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integrated  with  respect  to  Z  from  the  bottom  to  the  top  of  the  atmosphere.  Each  term 
in  the  equations  is  the  fmite  sum  of  separate  integrals.  Only  the  integrals  of  overlapping 
basis  functions  are  nonzero.  The  resultant  equations,  listed  in  Appendix  B,  are  matrix 
equations.  For  an  n-layer  model,  the  vectors  A,,  A2,  D,,  Dj,  H,,  Hj,  Q„  Qj,  T,,  Tj,  U,, 
Uj,  V,  and  Vj,  contain  n  +  2  components.  The  vectors  W,  and  Wj  contain  n+  1  com- 
ponents. The  matrices  M,  K,  and  O,  defmed  below,  are  (n  +  2)  x  (n  +  2)  matrices.  The 
matrix  P,  defmed  below,  is  an  (n+  1)  x  (n  +  2)  matrix.  The  mass  matrix  M  for  this 
model  is  defmed  bv 


M,j  =  f    (^j(Z)0.(Z)dZ  for  |i  -  jl  ^  1.  (2.46) 

The  matrix  K  is  defmed  for  terms  multiplied  by  u, 


1+1 
Kij(u)  =  \    u'^j%j(Z)0k(Z)(/),(Z)dZ  for  li  -  jl  <  1.  (2.47) 

k=i-l 

The  matrix  P  is  defined  for  terms  multiplied  bv  — --  W,  or  by  -:^:r-  W, 

'    dZ  cZ 


i+l 


P,3(x)  =    )^        -^  .A,(Z)(/).(Z)dZ  for  |i  -  jl  <  1, 

k=i-l      "^ 

where  x  is  u  or  T. 
The  matrix  0  is  defmed  for  terms  multiplied  by  -rr-  V, 


1+1 
<I>,j(u)  =    >    u*^       -^  cl>,iZ)HZ)^Z  for  |i  -  jl  ^  1. 

k=i-l        ^ 


(2.48) 


(2.49) 


The  staggered  basis  functions  present  two  general  problems  for  evaluating  the 
elements  of  the  four  matrices.  First,  for  an  n-layer  model,  portions  of  basis  functions 
4>i{Z)  and  ^n.2(Z)  are  defmed  in  the  model  atmosphere  but  the  physical  meaning  of 
contributions  from  those  terms  is  unclear.    The  contributions  are  included  in  the  first 
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two  rows  and  the  last  row  of  each  matrix.  Second,  only  portions  of  basis  functions 
02(Z)  and  ^n.i(Z)  are  defined  in  the  model  atmosphere.  To  describe  the  incomplete  sides 
of  both  basis  functions  an  assumption  must  be  made  about  the  value  oi'<p2  at  the  surface 
and  (/)„^i  at  the  top  of  the  atmosphere. 

Assumptions  are  made  and  procedures  are  developed  in  an  attempt  to  resolve 
these  problems.  In  this  model,  the  mean  state  variables,  ubar  and  Tbar,  are  defined  only 
at  the  n  staggered  levels.  However,  ubar  and  Tbar  values  defined  at  the  nodal  points 
of  </>i(Z)  and  0n^2(Z)  are  important  in  the  Galerkin  formulation  of  the  dubar/dZ  and 
dTba.T/8Z  terms.  In  these  experiments,  the  values  of  ubar  and  Tbar  are  defined  at  the 
surface  and  top  of  the  atmosphere.  Jordan  (1985)  did  not  define  them  at  the  nodal 
points  of  4>iiZ)  and  4>„,2{^)-  One  of  Shapiro's  (1987)  major  modifications  of  these  ex- 
periments was  to  define  ubar  and  Tbar  at  the  nodal  points  of  0,(Z)  and  ^^^^(Z).  For 
constant  shear  with  height,  ubar  and  Tbar  are  defined  at  the  boundaries  such  that  the 
shear  in  the  two  half  layers  at  the  boundaries  is  the  same  as  the  shear  in  the  other  layers. 
To  evaluate  the  staggered  basis  functions  defined  in  the  layers  between  the  surface  and 
Z2,  and  Zn^i  and  the  top  of  the  atmosphere,  it  is  assumed  that  the  value  at  the  boundaries 
of  those  basis  functions  is  one-half  Thus,  three-fourths  of  the  basis  functions 
02(Z)  and  </>n-i(Z)  are  defined  in  the  model  atmosphere. 

The  equations  for  the  general  elements  of  the  four  matrices  are  evaluated  by 
substituting  into  equations  (2.46)  through  (2.49)  the  formulas  for  (/>,.i(Z),  0,(Z).  0,_](Z), 
t/'i_i(Z),  </',(Z),  i/'j_i(Z),  and  i/'j_2(Z),  in  terms  of  the  local  coordinate  (^  =  Z  —  Z,.  The 
equations  for  these  basis  functions  defined  for  levels  1,  2,  i,  and  n+  1  are  listed  in  Ap- 
pendix C.  The  matrices  were  evaluated  by  integrating  numerically  using  2  point 
Gaussian  Quadratures  as  in  Shapiro  (1987). 

The  vorticity,  divergence  and  thermodynamic  equations,  written  in  matrix  and 
vector  form,  are 

M  -^  =  M(  -fD,  -  /?V,)  -  mK(u)A2,  (2.50) 

M  -^  =  M(  -fD2  -  J?V2)  +  A/K(u)A„  (2.51) 

M  -^  =  M(fA,  -  /?U,  +  f?H,)  -  /^K(u)D2  -  mP(u)W2,  (2.52) 
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M  -^  =  M(rA2  -  j5U2  +  a^'h^)  +  mK(u)D,  -  /.P(u)W„  (2.53) 

M  — -^  =  -  ;iK(u)T2  +  -^  <I>(u)V,  -  P(T)Wi  +  MQ,  and  (2.54) 

M  -^  =  /.K(u)T,  +  -|-  <D(u)V2  -  P(T)W2  +  MQ^.  (2.55) 

Equations  (2.50)  through  (2.55)  are  simplified  by  multiplying  each  equation  by 
M-'  and  applying  the  Robert  filter.  Actually,  one  should  not  compute  the  inverse  of 
M.  Instead,  at  t  =  0  one  should  obtain  the  LU  factorization  of  M.  Thus,  at  each  time 
step  one  only  needs  to  foru^ard  and  back  solve  a  triangular  system.  The  matrices 
M-'K,  M-'P  and  M-'O  are  constants.  They  are  constructed  in  the  mitialization  sub- 
routine and  stored  for  use  in  the  forecast  subroutine.  The  matrices  are  multiplied  by  the 
appropriate  vectors  with  values  for  time  level  nAt.  The  resultant  forecast  equations  are 
vector  equations  and  the  forecast  value  for  the  i-th  vertical  level  is  the  sum  of  values  in 
the  i-th  location  of  each  vector  equation.  The  prognostic  equations  for  the  vorticity, 
divergence  and  potential  temperature  vectors  are 

A^,+,)  =  A,(,_i)  +  2At(  -fD,  -  /?V,  -  MM-'K(u)A2)(n)'  (2-56) 

A2(n+1)  =  ^2in-l)  +  2At(  -fDj  -  /^V^  +  A^M-^K(u)AO(n),  (2-57) 

L^l(n+l)  =  ^l{n-l) 

+  2At(fA,  -  i5Ui  +  M^H,  -  mM"^K(u)D2  -  ;xM"^P(u)W2)(n),  (2.58) 

^2(n+l)  =^  ^2(n-l) 

+  2At(fA2  -  PV2  +  /^^Hj  -  //M-'K(u)Di  +  ;/M"'P(u)Wi)(„),  (2.59) 

*  l(n+1)  =  '^Un-l) 

+  2At(  -  mM-'K(u)T2  +  -^  M-'0(u)Vi  -  M-^P(T)W,  -f-  Q^\^^  and  (2.60) 
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+  2Al(/.M-'K(u)T,  +  -^  M-'<l>(u)V2  -  M-'P(T)W2  +  Q2)(n),  (2.61) 

where  the  subscripts  (n+  1),  (n)  and  (n-1)  refer  to  the  values  of  the  vectors  at  time  step 
(n+l)At,  nAt  and  (n-l)At,  respectively.  The  surface  geopotential  and  the  diagnostic 
variables  are  calculated  using  the  corresponding  equations  in  model  FDM-A  (see  Ap- 
pendix A). 

2.     FEM-B 

The  FEM-B  model  defines  vertical  velocity,  potential  temperature,  mean  state 
potential  temperature  and  diabaiic  heating  at  the  unstaggered  levels  in  terms  of  the  basis 
functions  )/')(Z).  The  other  variables  are  defined  at  the  staggered  levels  in  terms  of  the 
basis  functions  (/'j(Z).  The  basis  functions  are  the  same  as  defined  for  the  FEM-A  model, 
shown  in  Fig.  2. 

The  finite  element  approximations  for  the  vorticity,  divergence  and 
thermodynamic  equations  are  derived  by  substituting  the  expansion  for  each  dependent 
variable  into  equations  (2.25)  through  (2.30).  The  vorticity  and  divergence  equations 
are  multiplied  by  4),{Z)  and  integrated  with  respect  to  Z  from  the  bottom  to  the  top  oC 
the  atmosphere.  The  resultant  Galerkin  formulation  of  the  vorticity  and  divergence 
equations  are  the  same  as  those  derived  for  model  FEM-A.  The  matrices  in  those 
equations,  M,  K  and  P  are  given  by  equations  (2.46)  through  (2.48)  as  defined  for 
FEM-A.  The  thermodynamic  equations  are  multiplied  by  ip,{Z)  because  potential  tem- 
perature is  defined  at  the  unstaggered  levels.  As  before,  the  equations  are  integrated 
through  the  depth  of  the  atmosphere.  The  resultant  equations  are  listed  in  Appendix 
D.  Four  additional  matrices  are  defined  for  the  two  thermodynamic  equations.  The 
mass  matrix  11   is 

n,  =  [    ^/Z).A.(Z)dZ  for  |i  -  jl  ^  1.  (2.62) 

The  matrix  F  is  defined  for  terms  multiplied  by  u  , 
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i+1 


r,(u)  =  y  u'pii^^{Z)4>^{Z)il^,{Z)dZ  for  |i  -  jl  <  1.  (2.63) 

k=i-l 

Terms  multiplied  by  — —  V,  give  rise  to  the  transpose  of  the  matrix  P,  defined  by  (2.48). 
The  matrix  *F  is  defined  for  terms  multiplied  by  -z:=-  W, 


1+1 


Z-k       diAk 
t'       -^  ^j(Z).A,(Z)dZ  for|i-j|<l. 


(2.64) 

k=i-l       '^ 

As  discussed  in  the  FEM-A  model  description,  the  staggered  fmite  elements 
present  problems  for  evaluating  the  elements  of  the  matrices.  In  this  model,  ubar  is  de- 
fined at  the  surface,  the  top  of  the  atmosphere  and  at  the  n  staggered  levels.  The  mean 
state  temperature,  Tbar,  is  defined  at  the  unstaggered  levels  so  special  definitions  for  it 
are  not  needed.  Jordan  (1985)  did  not  include  the  contributions  from  the  perturbation 
quantities  defined  at  the  nodal  points  of  0,(Z)  and  (t>n+2(Z)-  They  were  included  in  this 
model  as  part  of  the  Shapiro  (1987)  modifications.  The  staggered  basis  functions, 
(t>,iZ),  are  evaluated  at  the  boundaries  using  the  assumptions  discussed  in  the  previous 
section. 

The  elements  of  matrices  11,  F,  and  T  are  evaluated  by  substituting  formulas 
for  4>,^i{Z),  0,_i(Z),  </),(Z),  </),_i(Z),  »/',^i(Z),  i//,(Z),  and  </'i-i(Z)  defined  in  terms  of  the  local 
coordinate  (*  =  Z  -  Z„  into  equations  (2.62),  (2.63)  and  (2.64).  Formulas  for  these  basis 
functions  are  hsted  in  Appendix  C.  As  in  FEiM-A,  the  matrices  are  evaluated  by  inte- 
grating numerically  using  2  point  Gaussian  Quadratures. 

The  forecast  matrix  equations  for  vorticity,  divergence  and  temperature  are  as 
described  for  model  FEM-A.  The  final  form  of  the  vorticity  and  divergence  vector 
equations  are  the  same  as  for  model  FEM-A,  equations  (2.56)  through  (2.59).  The 
thermodynamic  vector  equations  are 

Ti{n+1)  =  '^Un-l) 

+  2At(  -nU-'TiWi  +  -^  n-^P''(u)Vi  -  n-^4'(T)W,  +  Q,)(n).  (2.65) 
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+  2At(Mn-'r(u)T,  +  -^  n-'P^(u)V2  -  n-'4'(T)W2  +  Q^)^,),  (2.66) 

In  this  model,  the  vectors  Ai,  Aj,  D„  Dj,  H,,  H2,  Ui,  U2,  V,,  V2,  W,  and  Wj, 
contain  n  +  2  components.  The  vectors  Q,,  Q2,  T,  and  Tj  contain  n+1  components. 
Shapiro  (1987)  defined  vectors  W,  and  Wj  with  n+  1  components;  however,  since  W,  and 
Wj  are  derived  using  vectors  D,  and  Dj  (which  have  n+2  components)  in  equations 
(2.39)  and  (2.40),  W,  and  Wj  have  been  increased  to  n  +  2  components.  The  matrices 
n,  r,  and  T,  are  (n  +  1)  x  (n  +  1)  matrices  and  the  matrix  P'''  is  an  (n  +  2)  x  (n  +  1) 
matrix.  The  surface  geopotential  and  the  diagnostic  variables  are  calculated  using  the 
corresponding  equations  in  model  FDM-B. 
3.     FEM-C 

The  unstaggered  FEM  model  is  the  simplest  of  the  three  FEM  models.  Each 
of  the  dependent  variables  is  expanded  in  terms  of  the  basis  functions  0,(Z).  The  ex- 
pansion for  a  typical  term  is 

n+l 

A,(Z,t)=\  Ai(t)</)j(Z).  (2.67) 

j=i 

The  basis  functions  for  this  model  are  depicted  in  Fig.  3. 

The  finite  element  approximations  for  the  vorticity,  divergence  and 
thermodynamic  equations  are  derived  by  substituting  the  expansion  for  each  dependent 
variable  into  equations  (2.25)  through  (2.30).  Each  equation  is  multiplied  by  </),(Z)  and 
integrated  with  respect  to  Z  from  the  bottom  to  the  top  of  the  atmosphere.  Each  term 
in  the  equations  is  the  finite  sum  of  separate  integrals.  Only  the  integrals  of  overlapping 
basis  functions  are  nonzero.  The  resultant  equations,  listed  in  Appendix  E,  are  matrix 
equations.  The  matrices  M  and  K  are  given  by  equations  (2.46)  and  (2.47)  as  defined  for 
FEM-A.     The  matrix  E  is  defined  for  terms  multiplied  by  "T^  V,  by  —r^W,  or  by 

oZ       ' 


17 


Fig.  3.      Basis  functions  for  grid  C. 


E,j(x)  = 


i+1 

I 

k=i-l 


rz 


x*^ 


dZ 


0.(Z)0,(Z)dZ         forii-j|<l, 


Zo 


(2.68) 


where  x  is  u  or  T. 

Note  that  this  is  the  same  as  equation  (2.49)  which  is  only  defined  for  x  =  u  and  is  an 
(n  +  2)  X  (n  +  2)  matrix.  It  has  been  given  a  new  name  because  it  combines  equations 
(2.48)  and  (2.49)  and  is  now  an  (n  +  1)  x  (n  +  1)  matrix.  The  matrices  M,  K  and  E  are 
evaluated  using  the  same  procedures  as  models  FEM-A  and  FEM-B  except  only  the 
formulas  for  (f)^^)  are  needed  from  Appendix  C. 

The  vorticity  and  divergence  equations,  written  in  matrix  and  vector  form,  are 
given  in  equations  (2.50)  through  (2.53)  as  defmed  for  model  FEM-A.  The 
thermodynamic  matrix  equations  are 


1  T^ 

M  ^^  =  -  /iK(u)T2  +  -5-  E(u)V,  -  E(T)\V,  +  MQ,, 


dt 


R 


(2.69) 
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M  ^  =  A^K(u)Ti  +  ^  E(u)V2  -  E(T)W2  +  MQj.  (2.70) 

These  forecast  matrix  equations  are  simplified  as  described  for  model  FEM-A. 
The  final  form  of  the  vorticity  and  divergence  vector  equations  are  the  same  as  for  model 
FEM-A,  equations  (2.56)  through  (2.59).   The  thermodynamic  vector  equations  are 

Ti(n+1)  —  T^l(n-l) 

+  2At(  -  /^M"'K(u)T2  +  -|-  M"'E(u)V,  -  M"'E(T)W,  +  Qi)(„)  and  (2.71) 

'^2(n+l)  =  "^2(0-1) 

+  2At(A.M-'K(u)T,  +  ^  M-'E(u)V2  -  M-'E(T)W2  +  (^-,\,^_  (2.72) 

In  this  model,  the  vectors  A„  Aj,  D,,  D^,  H,,  H2,  Q,,  Q2,  T,,  T^,  Ui,  Uj,  Vi,  V^, 
W,  and  Wj  contain  n+  1  components.  The  matrices  M,  K  and  E  are  (n  +  1)  x  (n  4-  1) 
matrices.  The  surface  geopotential  and  the  diagnostic  variables  are  calculated  using,  the 
corresponding  equations  in  model  FDM-C. 
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III.     EXPERIMENTS  AND  RESULTS 

Several  experiments  are  performed  for  baroclinic  instability  with  all  or  some  of  the 
models.  The  first  results  use  the  linear  wind  profile  which  was  used  by  Shapiro  (1987). 
Most  of  the  experiments  use  a  hyperbolic  tangent  vertical  wind  profile  with  three  vari- 
ations. The  basic  experiments  use  an  initial  v-field  which  is  independent  of  Z  with  no 
temperature  perturbation.  These  experiments  are  integrated  for  96  hours  to  determine 
the  exponential  gro"wth  rate.  The  basic  experiments  use  two-,  four-,  six-  and  sixty-layer 
models.  The  next  to  last  experiment  uses  the  finite  difference  method  of  calculating  the 
thermodynamic  variables  in  the  finite  element  model  for  grid  A.  The  last  experiment 
uses  an  initial  v-field  which  tilts  backward  (toward  negative  x)  with  height.  The  temper- 
ature is  chosen  to  satisfy  the  thermal  wind  relation.  These  profiles  are  integrated  out  to 
72  hours  to  investigate  the  enhanced  initial  gro\\th  due  to  the  tilt  (Farrell,  1985).  The 
last  two  experiments  use  only  the  sixty-layer  models  except  for  the  thermodynamic  ex- 
periment which  includes  six-layer  models.  The  anahtical  solution  of  each  experiment 
does  not  exist  except  for  the  Eady  (1949)  solution  for  the  linear  shear  solutions. 

As  in  Jordan  (1985)  and  Shapiro  (1987),  the  sixty-layer  model  results  are  intercom- 
pared  to  determine  if  the  models  are  converging  to  the  same  solution.  The  standard  of 
comparison  for  the  two-,  four-  and  six-layer  models  is  the  consensus  sixty-layer  solution. 
The  v-component  and  temperature  profiles  are  examined  in  each  experiment.  Expo- 
nential grovnh  rates  are  calculated  to  help  compare  the  models. 

Several  parameters  are  defined  identically  in  each  experiment.  The  vertical  coordi- 
nate, Z,  is  defined  between  zero  and  one  (1000-368mb)  and  the  vertical  levels  are  equally 
spaced.  The  Coriolis  parameter  is  defined  at  45  degrees  latitude.  The  mean  state  po- 
tential temperature  increases  with  height  from  its  surface  value  of  310.0  degrees  Kelvin. 
The  vertical  shear  in  the  ubar  field  is  defined  in  each  experiment  as  a  function  of  Z 
multipUed  by  STRGTH,  where  STRGTH  defines  the  strength  of  the  wind  at  the  top  of 
the  atmosphere  in  m/s.  In  all  experiments,  STRGTH  is  defmed  as  40  m/s.  There  is  no 
diabatic  heat  source  and  no  mountain  topography.  Unless  otherwise  indicated,  all 
forecast  experiments  are  terminated  at  96  hours. 
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A.     LINEAR  VERTICAL  SHEAR  IN  THE  UBAR  FIELD 

The  wind  profile  in  these  experiments  is  a  linear  function  of  Z,  with  ubar  = 
{STRGTH)Z.   Waves  are  generated  in  each  model  using  an  initial  perturbation  ofv'  = 
5.0  m.  s  in  the  cosine  term.    Wavelengths  (in  x)  of  4000,  3000  and  2000  km  are  used  for 
each  model  for  two,  six  and  sixty  levels. 
1.    Sixty- Layer  Models 

For  an  x-wavelength  of  4000  km,  the  sixty-layer  FEM-C  and  FDM-C  models 
both  converge  to  the  same  solution  as  the  FDM-A  and  FDM-B  models  from  Shapiro 
(1987)  for  the  v-component  amplitude  and  phase  (Fig.  4)  and  the  temperature  amplitude 
(Fig.  5).  All  figures  will  be  found  at  the  end  of  the  chapter.  As  in  Shapiro  (1987),  the 
phase  is  defined  between  zero  and  360  degrees  and  a  discontinuity  occurs  in  the  phase 
profile  if  the  phase  passes  through  zero  degrees.  FDM-A,  FDM-B,  FDM-C  and  FEM-C 
represent  the  same  physical  solution,  which  is  called  the  consensus  solution.  These  sol- 
utions all  have  the  same  patterns  which  were  derived  by  Eady  (1949)  for  the  quasi- 
geostrophic  model.  The  v-field  has  a  maximum  amplitude  on  the  boundaries  and  tilts 
backward  with  height.  The  FEM-B  amplitudes  in  Fig.  4  and  Fig.  5  are  less  than  the 
consensus  which  indicates  a  smaller  growth  rate.  The  FEM-B  temperature  field  (Fig. 
5)  also  has  an  oscillation  near  the  surface  similar  to  what  was  found  by  Shapiro  even 
though  an  additional  element  is  added  to  the  vertical  velocity.  The  FEM-A  scheme  has 
a  slightly  larger  gro\^th  rate  than  the  consensus.  This  suggests  that  the  FEM-A  and 
FEM-B  models  need  further  improvement  to  take  full  advantage  of  their  possibilities. 

For  an  x-wavelength  of  3000  km,  the  sixty-layer  profiles  for  the  v-component 
and  temperature  show  the  same  results  as  the  4000  km  case  except  that  the  magnitude 
of  the  amphtudes  are  smaller  and  the  profiles  appear  to  contain  some  oscillation  (Fig. 
6).  The  reduced  amplitude  agrees  with  the  Eady  (1949)  theor}"  which  gives  smaller 
grovnh  rates  for  this  wavelength. 

At  a  smaller  wavelength  of  2000  km,  the  v-component  and  temperature  ampli- 
tudes show  verv'  little  growth  at  all  but  have  several  oscillations  throughout  the  atmos- 
phere (Figs.  7-8).  The  finite  difference  models  and  FEM-C  converge  to  the  same 
consensus  solution.  FEM-A  is  very  close  to  the  consensus  profile,  but  has  a  shghtly 
smaller  amplitude  in  the  middle  layers  of  the  temperature  profile  and  a  slightly  larger 
amplitude  in  the  v-component  profile.  FEM-B  has  a  smaller  amplitude  than  the  con- 
sensus profile  for  the  v-component  and  oscillates  between  having  a  smaller  amplitude  to 
having  a  larger  amplitude  in  the  temperature  profile.  The  phase  profiles  show  that  all 
the  models  except  FEM-A  and  FEM-B  converge  to  the  same  phase  solution  (Fig.  7). 
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FEM-A  is  ven"  close  to  the  consensus  solution  but  changes  from  being  slightly  greater 
in  the  lower  layers  to  being  slightly  less  in  the  upper  layers.  FEM-B  is  smaller  than  the 
consensus  in  the  lower  layers  and  larger  in  the  upper  layers.  The  lack  of  growth  in  this 
case  agrees  with  the  Eady  theor\^  which  gives  only  neutral  modes  for  this  scale.  Farrell 
(1985)  has  shown  that  a  given  initial  state  can  funnel  energy  into  the  neutral  modes  of 
the  problem.  In  this  case,  the  most  likely  neutral  waves  which  appear  in  the  solutions 
are  the  two  quasi-geostrophic  solutions  and  internal  gravity  waves  with  small  vertical 
scales. 

2.     SL\- Layer  Models 

For  the  4000  km  wavelength,  all  the  models  except  FEM-B  and  FEM-A  con- 
verge to  the  consensus  for  the  v-component  amplitude  (Fig.  9)  and  all  the  models  except 
FEM-B  are  converging  to  the  consensus  solution  for  the  temperature  amplitude  (Fig. 
10).  FEM-B  has  the  same  shape  as  the  consensus  profile  but  has  a  smaller  amplitude 
for  temperature  and  the  v-component.  FEM-A  has  a  much  smaller  v-component  than 
the  consensus.  The  v-component  phase  profiles  show  that  all  the  models  except  FEM-A 
and  FEM-B  converge  to  the  consensus  phase  profile  (Fig.  9).  The  six-layer  FDM-B 
model  is  closest  to  the  consensus. 

At  3000  km,  the  difference  among  the  six-layer  models  and  between  the  six-layer 
models  and  the  consensus  profile  for  v-component  and  temperature  amplitudes  becomes 
more  pronounced  (Figs.  11-12).  In  general,  the  finite  difference  models  are  closer  than 
the  finite  element  models  to  the  consensus  profile  for  schemes  B  and  C.  The  opposite 
is  true  for  grid  A  except  in  the  lower  layers  of  the  temperature  amplitude  profile. 
FEM-B  has  a  much  smaller  amplitude  than  the  consensus.  The  phase  profiles  show  the 
models  are  still  converging  to  the  consensus  phase  profile  except  for  FEM-B  (Fig.  11). 

At  2000  km,  the  models  appear  to  be  a  little  chaotic  (Figs.  13-14).  However, 
all  the  models  depict  the  overall  shape  of  the  consensus  profile  except  for  the  oscillations 
in  the  temperature  amplitude  profiles.  FDM-B  appears  to  represent  the  consensus  better 
than  the  other  models  although  FEM-B  gives  the  same  surface  temperature  prediction 
as  the  consensus  profile.  In  the  v-component  amplitude  profiles,  the  finite  element 
models  are  closer  then  the  finite  difference  models  to  the  consensus  profile  for  grids  A 
and  C  but  not  B.  However,  the  surface  prediction  for  FDM-C  is  closer  than  the  pre- 
diction for  FEM-C.  FDM-B  is  closest  to  the  consensus  profile  for  the  v-component 
amplitude.  The  phase  profiles  show  that  the  models  for  scheme  A  are  closer  to  the 
consensus  phase  profile  than  are  the  models  for  schemes  B  and  C  (Fig.  13).  Clearly,  the 
small  oscillations  can  not  be  reproduced  in  any  six-layer  scheme. 
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3.     Two-Layer  Models 

The  curves  for  the  two-layer  models  are  not  shown  to  save  space.  For  an  x- 
wavelength  of  4000  km,  FDM-C  represents  the  shape  of  the  consensus  profile  of  tem- 
perature and  v-component  better  than  the  other  models,  but  it  has  the  benefit  of  an 
extra  prediction  level  for  these  variables.  FEM-C  gives  the  best  surface  prediction  for 
both  the  temperature  and  v-component.  In  the  v-component  amplitude  profiles,  the  fi- 
nite difierence  models  are  closer  to  the  consensus  amplitude  than  are  the  finite  element 
models  for  schemes  A  and  B.  However,  these  two  schemes  do  not  have  enough  resol- 
ution to  depict  the  shape  of  the  consensus  profile.  The  phase  profiles  for  the  v- 
component  show  general  agreement  among  the  models  in  representing  the  shape  of  the 
consensus. 

At  3000  km,  FDM-C  again  represents  the  shape  of  the  consensus  amplitude 
profile  for  temperature  and  v-component  well.  FEM-A  is  also  close  to  the  magnitude 
of  the  consensus  profile.   The  other  models  give  poorer  results. 

At  2000  km,  the  magnitude  of  FDM-C  overshadows  all  the  other  models. 
FEM-C  and  FDM-A  also  show  a  tendency  for  growlh.  In  the  v-component  phase 
profiles,  all  the  A  and  B  models  have  a  negative  tilt  instead  of  the  positive  tilt  of  the 
consensus  profile.  FDM-C  has  a  positive  tilt  in  the  lower  atmosphere  and  FEM-C  has 
a  positive  tilt  in  the  upper  atmosphere. 

B.     NONLINEAR  VERTICAL  SHEAR  IN  THE  UBAR  FIELD 

The  nonlinear  vertical  shear  is  created  by  using  a  hyperbolic  tangent  function  of  Z 
similar  to  the  one  used  by  Satyamurty  et  al  (1982).  In  this  experiment,  ubar  = 
(STRGTH)  tanh(aZ  +  K)  where  a  and  k  are  given  difierent  values  to  control  the  magni- 
tude and  location  of  the  shear.  Fig.  15  illustrates  the  three  ubar  profiles  used  in  this 
experiment  which  are  referred  to  as  Type  I,  Type  II  and  Type  III.  Waves  are  generated 
in  each  model  using  an  initial  perturbation  of  v'  =  5.0  m/s  in  the  cosine  term.  The  x- 
wavelengths  of  4000,  3000  and  2000  km  are  used  for  each  model. 
1.     Hyperbolic  Shear  Type  I 

For  this  part  of  the  experiment,  a  =  2.0  and  k  =  0.0.  This  increases  the  mag- 
nitude of  the  low  level  ubar  shear  and  forces  the  largest  gradient  to  occur  near  the  sur- 
face. 

a.     Sixty-Layer  Models 

Table  1  contains  the  exponential  growth  rates  for  all  the  sixty-layer  models. 
For  an  x-wavelength  of  4000  km,  the  growth  rates  for  Type  I  indicate  that  the  sixty-layer 
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finite  diflerence  models  and  FEM-C  all  converge  to  the  same  growth  rate  after  96  hours. 
This  is  similar  to  the  results  from  the  linear  wind  profile. 


Table  1 

Exponential  Grovnh  Rates  (  x  10-^) 

Hyperbolic  Shear 

60-layer  Models 

FDM-A 

FDM-B 

FDM-C 

FEM-A 

FEM-B 

FEM-C 

Type  I 

4000  km 
3000  km 
2000  km 

1.45 
1.53 
1.33 

1.45 
1.53 
1.33 

1.45 
1.53 
1.33 

1.46 
1.54 
1.35 

1.43 
1.48 
1.28 

1.45 
1.53 
1.32 

Type  II 

4000  km 
3000  km 
2000  km 

1.82 
2.14 
2.46 

1.81 
2.14 
2.46 

1.82 
2.14 

2.46 

1.82 

2.14 
2.46 

1.78 
2.09 
2.40 

1.82 

2.14 
2.47 

Type  III 

4000  km 
3000  km 
2000  km 

2.89 
3.21 
3.02 

2.89 

3.22 
3.03 

2.89 

3.22 
3.03 

2.89 
3.21 
3.01 

2.83 
3.14 
2.90 

2.89 

3.22 
3.04 

The  growth  rate  for  FEM-A  is  slightly  larger  and  the  growth  rate  for  FEM-B  is  slightly 
smaller.  FDM-A,  FDM-B,  FDM-C  and  FEM-C  all  converge  to  the  same  solution  for 
the  v-component  and  temperature  amplitudes  and  represent  the  consensus  solution.  The 
FEM-A  model  has  a  slightly  higher  amplitude  than  the  consensus  profile  and  the  shape 
is  the  same  as  the  consensus  profile.  The  FEM-B  model  has  a  lower  temperature  am- 
plitude than  the  consensus  profile,  but  the  shape  is  the  same  except  for  an  oscillation  in 
the  lower  layers  similar  to  the  linear  shear  case.  The  difference  between  FEM-B  and  the 
consensus  profile  is  greater  in  the  lower  layers  than  in  the  upper  layers.  For  the  v- 
component  amplitude,  the  amplitude  for  FEM-B  is  again  lower  than  the  consensus 
profile,  but  there  is  no  oscillation  in  the  lower  layers  of  the  atmosphere.  The  profiles  for 
the  v-component  phase  are  the  same  for  all  six  models.  The  profiles  are  not  shown  be- 
cause the  consensus  profile  is  included  with  the  lower  resolution  curves. 
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For  an  x-wavelength  of  3000  km,  the  growth  rates  show  the  same  the  results 
as  4000  km  except  the  growth  rates  are  greater.  The  v-component  and  temperature 
profiles  also  show  the  same  results  as  4000  km.  Again,  the  amplitudes  are  larger  for  the 
smaller  wavelength.  The  v-component  phase  profiles  show  that  FDM-A,  FDM-B, 
FDM-C  and  FEM-C  are  identical,  but  FEM-A  and  FEM-B  have  a  slightly  larger  phase. 

For  an  x-wavelength  of  2000  km,  the  growth  rates  show^  the  same  results 
as  4000  km  except  the  growth  rates  are  smaller.  The  v-component  and  temperature 
profiles  also  show  the  same  results  as  4000  km,  but  the  amplitudes  are  smaller.  The 
ampUtude  profiles  also  appear  to  contain  some  oscillation.  The  v-component  phase 
profiles  are  identical  for  all  the  models  except  FEM-A  which  has  a  slightly  larger  phase. 
This  profile  clearly  has  the  instabiUty  shifted  to  smaller  wavelengths  than  the  linear 
profile  due  to  the  vertical  wind  shear  being  larger  near  the  surface. 
b.     Six-layer  Models 

The  profiles  for  the  six-  and  two-layer  models  are  not  shown  to  save  space. 
Table  2  contains  the  exponential  growth  rates  for  six-,  four-  and  two-layer  models  for 
Type  I.  For  an  x-wavelength  of  4000  km,  the  growth  rates  show^  that  all  the  models 
except  FEM-B  are  converging  to  the  growth  rate  of  the  consensus  profile.  The  v- 
component  amplitude  profiles  support  the  convergence  of  the  models  to  the  consensus 
profile.  The  finite  difference  models  are  closer  to  the  consensus  solution  than  are  the 
finite  element  models  for  grids  A  and  B,  but  FEM-C  is  closer  than  FDM-C.  All  models 
represent  the  low  level  shear  well.  Although  the  growth  rates  show  that  both  FDM-A 
and  FEM-C  are  closest  to  the  consensus,  the  amplitude  profiles  show*  that  FEM-C  is 
closest  to  the  consensus  profile.  The  temperature  amplitude  also  shows  the  convergence 
of  the  models  to  the  consensus  solution  except  for  FEM-B.  In  this  case,  the  surface 
prediction  for  both  FEM-B  and  FDM-B  is  considerably  less  than  the  consensus.  The 
v-component  phase  profiles  show  a  moderate  amount  of  negative  tilt  in  the  lower  half 
of  the  atmosphere  corresponding  to  the  low  level  shear  region.  All  the  models  have  the 
same  phase  profile  shape,  but  scheme  C  is  closer  to  the  consensus  solution  than  are 
schemes  A  and  B.  This  difference  is  due  to  a  difference  in  phase  speed  between  the 
schemes. 
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Table  2 

Exponential  Growth  Rates  (  x  10"-) 

Hyperbolic  Shear  Type  I 

Models 

Consensus 

FDM-A 

FDM-B 

FDM-C 

FEM-A 

FEM-B 

FEM-C 

6  Layers 

4000  km 
3000  km 
2000  km 

1.45 
1.53 
1.33 

1.44 
1.52 

1.37 

1.42 
1.48 
1.41 

1.49 
1.54 
1.46 

1.43 
1.49 

1.37 

1.36 
1.41 
1.52 

1.46 
1.55 

1.54 

4  Layers 

4000  km 
3000  km 
2000  km 

1.45 
1.53 
1.33 

1.35 
1.40 
1.50 

1.36 

1.48 
1.4S 

1.58 
1.66 

-0.20* 

1.32 
1.31 

1.15 

1.33 

1.46 
1.55 

1.51 
1.64 

1.43 

2  Layers 

4000  km 
3000  km 
2000  km 

1.45 
1.53 
1.33 

1.29 
1.51 
1.75 

1.03 

1.00 

0.24* 

1.49 
1.83 
2.66 

1.09 
1.13 
1.10 

1.07 
0.98 

-0.87* 

1.36 
1.51 

2.07 

*  These  c 

ases  do  not  c 

onverge  to  a  steady  growth  rate 

The  gro\«h  rates  for  an  x-wavelength  of  3000  km  give  results  similar  to 
those  for  4000  km  except  that  FDM-C  is  closer  than  FEM-C  to  the  consensus  solution. 
The  v-component  amplitude  and  phase  profiles  are  also  similar.  In  this  case,  however, 
FDM-A  is  closest  to  the  consensus  amplitude  profile  instead  of  FE.M-C,  but  FEM-B  is 
still  farthest  from  the  consensus.  The  phase  results  are  identical  to  the  results  for  the 
4000  km  wavelength  case. 

At  an  x-wavelength  of  2000  km,  a  few  more  changes  are  observed  from  the 
longer  wavelengths.  For  this  wavelength,  the  growth  rates  indicate  that  all  the  six-layer 
models  are  growing  faster  than  the  consensus  solution.  Although  the  growth  rates  for 
FDM-A  and  FEM-A  are  the  same,  the  difference  between  the  finite  difference  and  finite 
element  models  for  grids  B  and  C  is  much  greater  at  the  smaller  wavelength.  These  dif- 
ferences are  also  observed  in  the  v-component  and  temperature  amplitude  profiles.  Al- 
though the  grovnh  rates  for  FDM-A  and  FEM-A  are  the  same,  the  amplitude  profiles 
show  that  FEM-A  is  slightly  closer  to  the  v-component  consensus  solution,  but  FDM-A 


26 


is  slightly  closer  to  the  temperature  amplitude  consensus.    This  may  occur  because  the 
growth  rate  is  computed  at  the  first  layer  above  Z  =   0  and  the  difTerences  between 
FDM-A  and  FEM-A  occur  above  this  level. 
c.     Four-Layer  Models 

In  the  4000  km  wavelength  case,  the  growth  rates  show  that  the  four-layer 
models  do  not  converge  to  the  consensus  solution  as  well  as  the  six-layer  4000  km 
models.  Otherwise,  the  results  are  basically  the  same.  The  finite  difTerence  models  are 
closer  than  the  finite  element  models  for  grids  A  and  B,  but  not  grid  C.  However,  the 
v-component  amplitude  profiles  indicate  that  FEM-C  is  closest  only  at  the  second  layer, 
where  the  growth  rate  is  computed,  and  the  top  layer  (Fig.  16).  The  amplitude  profiles 
also  show  the  lesser  convergence  of  the  four-layer  models  to  the  consensus  solution. 
The  v-component  phase  profiles  give  the  same  results  as  the  six-layer  models,  but  with 
a  slightly  greater  difference  between  the  four-layer  models  and  the  consensus. 

For  the  3000  km  wavelength,  the  growth  rates  for  the  four-layer  models 
again  indicate  results  similar  to  the  six-layer  models.  The  finite  difference  models  are 
closer  to  the  consensus  solution  than  are  the  fmite  element  models,  but  FDM-B  is  clos- 
est to  the  consensus  growth  rate  instead  of  FDM-A.  FEM-A  has  the  worst  growth  rate 
in  this  case.  Also,  the  four-layer  models  do  not  converge  to  the  consensus  solution  as 
well  as  do  the  six-layer  models  (Fig.  17).  The  v-component  amplitude  profiles  support 
the  growth  rate  results.  The  phase  results  for  the  four-layer  models  are  similar  to  the 
results  for  the  six-layer  models.  The  primary  difference  is  that  except  for  FEM-C  the 
four-layer  model  phase  profiles  are  not  as  close  to  the  consensus  phase  profile  as  are  the 
six-layer  models. 

For  the  2000  km  wavelength,  the  growth  rates  indicate  much  less  conver- 
gence of  the  four-layer  models  to  the  consensus  solution.  In  fact,  FDM-C  fails  to  con- 
verge to  a  steady  growth  rate  altogether.  For  this  case,  the  fmite  element  model  growth 
rates  are  closer  than  the  finite  difference  model  rates  to  the  consensus  solution  for 
schemes  A  and  C,  but  not  scheme  B.  All  the  models  except  FEM-A  and  FDM-C  grow 
much  faster  than  the  consensus.  FEM-C  is  closest  to  the  consensus  growth  rate. 
FEM-C  also  reflects  the  shape  of  the  v-component  amplitude  consensus  profile  better 
than  the  other  four-layer  models  (Fig.  18).  All  the  models  except  FDM-C  depict  the 
same  shape  as  the  phase  consensus  and  FEM-C  is  again  closest  to  the  consensus  sol- 
ution.  Note  that  the  shorter  waves  also  have  smaller  vertical  scales. 
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d.     Two- Layer  Models 

For  the  4000  km  wavelength,  the  growth  rates  indicate  that  the  two-layer 
finite  difierence  models  are  closer  to  the  consensus  solution  than  are  the  finite  element 
models  for  grids  A  and  C,  but  not  grid  B.  FEM-C  is  closest  to  the  consensus  amplitude 
profile  for  v-component  and  temperature.  Schemes  A  and  B  do  not  have  enough  resol- 
ution to  depict  the  low  level  shear  in  the  v-component  amplitude  profile.  The  grid  C 
models  attempt  to  reflect  the  low  level  shear,  but  the  surface  predictions  are  considerably 
less  than  the  consensus  solution.  All  six  models  reflect  the  negative  tilt  of  the  consensus 
phase  profile  even  though  they  have  only  two  levels.  FEM-C  is  closest  to  the  consensus 
phase  profile. 

The  growth  rates  for  the  3000  km  case  indicate  that  the  finite  difference 
models  are  closer  than  the  finite  element  models  to  the  consensus  growth  rate  for  grids 
A  and  B,  but  FEM-C  is  closer  than  FDM-C.  The  scheme  B  models  and  FEM-A  have 
much  slower  growih  rates  than  the  consensus  and  FDM-C  has  a  larger  growth  rate. 
FEM-C  is  closest  to  the  consensus  amplitude  profile  for  the  v-component,  but  FDM-A 
is  closest  in  the  middle  and  lower  levels  of  the  temperature  profile.  Although  grid  B  also 
predicts  the  temperature  at  the  surface,  its  prediction  is  much  poorer  than  the  consensus 
solution.  The  two-layer  A  and  B  models  do  not  have  enough  resolution  to  depict  the 
shape  of  consensus  v-component  phase  profile  but  they  still  reflect  the  negative  tilt. 

For  the  2000  km  wavelength,  the  growlh  rates  for  FDM-A,  FDM-C  and 
FEM-C  are  much  greater  than  the  consensus  growth  rate,  and  FEM-A  is  slightly  less 
than  and  closest  to  the  consensus.  FDM-B  and  FEM-B  grow  very  little  and  their  growth 
rates  do  not  converge  to  a  steady  gro\\th  rate.  The  finite  element  models  are  closer  than 
the  finite  difference  models  to  the  consensus  v-component  and  temperature  amplitude 
profiles  for  schemes  A  and  C,  but  not  for  scheme  B.  FEM-A  is  closest  in  magnitude  to 
the  consensus  solution. 

2.     Hyperbolic  Shear  Type  II 

For  this  part  of  the  experiment,  a  =  4.0  and  k  =  0.0.    This  increases  the  mag- 
nitude of  the  low  level  ubar  shear  but  keeps  the  largest  gradient  near  the  surface. 
a.     Sixty-Layer  Models 

For  an  x-wavelength  of  4000  km,  the  growth  rates  for  Type  II  indicate  that 
the  sixty-layer  A  and  C  models  all  converge  to  the  same  grov^th  rate  after  96  hours 
(Table  1).  The  growth  rate  for  FDM-B  is  nearly  the  same  but  slighly  smaller  and  the 
growth  rate  for  FEM-B  is  the  smallest.  FDM-A,  FDM-C,  FEM-A  and  FEM-C  all 
converge  to  the  same  solution  for  the  v-component  and  temperature  amplitudes  and 
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represent  the  consensus  solution.  The  FDM-B  model  has  a  slightly  smaller  amplitude 
than  the  consensus  profile,  but  the  shape  is  the  same  as  the  consensus  profile.  The 
FEM-B  model  has  a  lower  temperature  amplitude  than  the  consensus  profile  and  the 
shape  is  the  same  except  for  a  similar  oscillation  to  that  of  the  linear  shear  case  in  the 
lower  layers.  The  difference  between  FEM-B  and  the  consensus  profile  is  greater  in  the 
lower  layers  than  in  the  upper  layers.  For  the  v-component  amplitude,  the  amplitude 
for  FEM-B  is  again  lower  than  the  consensus  profile,  but  there  is  no  oscillation  in  the 
lower  layers  of  the  atmosphere.  The  magnitudes  of  the  growlh  rates  and  the  amplitude 
profiles  are  greater  than  for  the  Type  I  case.  The  profiles  for  the  v-component  phase 
are  the  same  for  all  six  models.  The  sixty-layer  profiles  are  not  shown  because  the  con- 
sensus profile  is  included  with  the  lower  resolution  curves. 

For  an  x-wavelength  of  3000  km,  the  growth  rates  for  all  the  sixty-layer 
models  are  identical  except  for  FEM-B  and  are  greater  than  the  4000  km  growth  rates. 
The  v-component  and  temperature  profiles  show  results  similar  to  the  4000  km  case,  but 
with  larger  amplitudes.  The  v-component  phase  profiles  show  that  all  the  models  have 
identical  phases  at  96  hours  except  for  FEM-A  which  has  a  slightly  larger  phase. 

For  an  x-wavelength  of  2000  km,  the  results  are  the  same  as  for  the  3000 
km  case  except  FEM-C  has  a  slightly  larger  growth  rate.  However,  the  v-component 
and  temperature  profiles  show  the  same  results  as  for  4000  km  except  the  amplitudes  are 
larger  for  2000  km.  The  v-component  phase  profiles  are  identical  for  all  the  models  ex- 
cept FEM-A  which  has  a  slightly  larger  phase. 
b.     Six- Layer  Models 

To  save  space,  the  profiles  for  the  six-  and  two-layer  models  are  not  shown. 
Table  3  contains  the  exponential  grovilh  rates  for  six-,  four-  and  two-layer  models  for 
Type  II.  For  an  x-wavelength  of  4000  km,  the  growth  rates  indicate  that  the  finite  ele- 
ment models  are  closer  to  the  consensus  growth  rate  than  are  the  finite  difference  models 
for  grids  B  and  C,  but  not  for  grid  A.  The  gro\\th  rates  for  grid  C  are  closer  than  the 
rates  for  grids  A  and  B  to  the  consensus.  FEM-C  is  closest  to  the  consensus  growth 
rate.  All  the  models  represent  the  shape  of  the  consensus  profile  well  except  for  the  grid 
B  temperature  prediction.  FEM-C  is  closest  to  the  consensus  amplitude  profile.  The 
v-component  phase  profiles  show  the  negative  tilt  in  the  lower  half  of  the  atmosphere 
corresponding  to  the  region  of  low  level  shear.  All  the  six-layer  models  have  the  same 
shape  as  the  consensus  phase  profile,  but  scheme  C  is  closer  to  the  consensus  than  are 
schemes  A  and  B. 
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Table  3 

Exponential  Growth  Rates  (  x  10"-) 

Hyperbolic  Shear  Type  II 

Models 

Consensus 

FDM-A 

FDM-B 

FDM-C 

FEM-A 

FEM-B 

FEM-C 

6  Layers 

4000  km 
3000  km 
2000  km 

1.82 

2.14 
2.46 

1.68 
1.95 

2.29 

1.63 
1.95 

2.43 

1.93 

2.26 
2.14 

1.66 
1.89 
2.05 

1.67 
1.99 

2.45 

1.87 

2.23 
2.54 

4  Layers 

4000  km 
3000  km 
2000  km 

1.82 

2.14 
2.46 

1.63 
1.97 

2.48 

1.36 
1.61 
1.82 

1.82 
1.89 

2.47 

1.51 
1.70 
1.98 

1.35 
1.59 
1.84 

1.78 
2.05 
2.25 

2  Layers 

4000  km 
3000  km 
2000  km 

1.82 

2.14 
2.46 

1.16 

1.38 

1.72 

0.57 

1.42 

0.25* 

1.80 
2.32 
3.16 

1.17 
1.23 
1.30 

0.49 

0.56* 

-0.34* 

1.58 
1.90 

2.57 

*  These  cases  do  not  converge  to  a  steady  groviah  rate 

For  the  3000  km  wavelength  case,  the  growth  rates  indicate  results  similar 
to  the  4000  km  case.  For  grids  A  and  C,  the  finite  difference  models  are  closer  to  the 
consensus  profile  than  are  the  finite  element  models  although  the  differences  are  greater 
in  the  low  level  shear  region  than  in  the  upper  levels.  FEM-B  is  closer  than  FDM-B  to 
the  consensus  amplitude  profile  except  at  the  midpoint  of  the  temperature  profile. 
FDM-C  is  closest  to  the  consensus  profile  except  in  the  lowest  levels  of  the  temperature 
profile  above  the  surface  where  FEM-B  is  closest.  The  phase  profiles  show  the  same 
results  as  the  4000  km  case.  For  the  2000  km  wavelength  case,  all  the  results  are  the 
same  as  the  results  for  4000  km. 
c.     Four- Layer  Models 

The  growth  rates  for  the  four-layer  4000  km  wavelength  models  indicate  the 
finite  difference  models  are  closer  than  the  finite  element  models  to  the  consensus  groAvth 
rate.    Schemes  A  and  B  have  much  smaller  growth  rates  than  both  scheme  C  and  the 
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consensus  solution.  FDM-C  has  the  same  growth  rate  as  the  consensus.  Overall,  the 
v-component  and  temperature  amplitude  profiles  show  that  the  finite  difference  models 
are  closer  to  the  consensus  profile  than  are  the  finite  element  models  (Figs.  19-20). 
However,  for  scheme  B  the  finite  element  model  is  closer  to  the  consensus  amplitude 
profile  for  the  temperature.  For  scheme  C,  the  finite  element  model  is  closer  to  the 
consensus  at  the  surface  and  at  Z  =  0.25  for  the  v-component  profile.  Otherwise, 
FDM-C  is  closest  to  the  consensus  profile.  The  v-component  phase  profiles  show  the 
negative  tilt  in  the  lower  half  of  the  atmosphere  corresponding  to  the  region  of  low  level 
shear  (Fig.  19).  All  the  six-layer  models  have  the  same  shape  as  the  consensus  phase 
profile,  but  scheme  C  is  closer  to  the  consensus  than  are  schemes  A  and  B. 

For  the  3000  km  wavelength,  the  growth  rates  for  the  four-layer  models  in- 
dicate results  similar  to  the  4000  km  case,  except  that  FEM-C  is  closer  than  FDM-C  to 
the  consensus  solution.  FEM-C  is  also  closest  to  the  growth  rate  for  the  consensus. 
Overall,  the  finite  difference  models  are  closer  than  the  fmite  element  models  to  the 
consensus  amplitude  profile  for  schemes  A  and  B,  but  not  for  scheme  C  in  the  v- 
component  profiles  (Fig.  21).  The  finite  element  models  are  closer  than  the  finite  dif- 
ference models  for  grids  B  and  C,  but  not  grid  A  for  the  temperature  profile  (Fig.  22). 
FEM-C  is  closest  to  the  consensus  except  in  the  middle  layers  of  the  temperature  profile 
where  FDM-A  is  closest.  For  the  phase  profiles,  the  results  for  this  case  are  the  same 
as  the  results  for  the  4000  and  3000  km  cases  (Fig.  21). 

Growth  rates  for  the  four-layer  models  in  the  2000  km  case  indicate  that  the 
finite  difference  models  are  closer  than  the  finite  element  models  to  the  consensus  sol- 
ution for  grids  A  and  C.  FDM-B  and  FEM-B  are  nearly  identical.  FDM-C  is  closest 
to  the  growth  rate  for  the  consensus  solution.  In  this  case,  the  finite  element  models  are 
closer  than  the  finite  difference  models  to  the  consensus  ampHtude  profile  for  scheme  B 
in  the  lower  half  of  the  atmosphere  and  for  scheme  C,  but  not  for  scheme  A  except  in 
the  upper  levels  of  the  v-component  amplitude  profile  (Fig.  23).  FEM-C  is  closest  to  the 
consensus  profile. 

d.     Two- Layer  Models 

For  the  4000  km  wavelength,  the  growth  rates  indicate  that  the  two-layer 
finite  difierence  models  are  closer  to  the  consensus  solution  than  are  the  finite  element 
models  for  grids  B  and  C,  but  not  for  grid  A.  FEM-C  is  slightly  closer  than  FDM-A  to 
the  consensus.  FDM-C  is  closest  to  the  growth  rate  for  the  consensus  solution,  but 
FEM-C  is  closest  overall  to  the  consensus  amplitude  profile  for  v-component  and  tem- 
perature.   Schemes  A  and  B  do  not  have  enough  resolution  to  depict  the  low^  level  shear 
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in  the  v-component  amplitude  profile.  The  grid  C  models  attempt  to  reflect  the  low  level 
shear,  but  the  surface  predictions  are  considerably  less  than  the  consensus  solution.  All 
six  models  refiect  the  negative  tilt  of  the  v-component  even  though  they  have  only  two 
levels  and  FEM-C  is  closest  to  the  consensus  phase  profile. 

For  the  3000  km  case,  the  gro\M:h  rates  show  that  the  finite  difference 
models  are  closer  than  the  finite  element  models  to  the  concensus  for  schemes  A  and 
C,  but  not  for  scheme  B.  FEM-B  appears  to  be  closer  than  FDM-B  to  the  consensus 
growth  rate,  but  FEM-B  does  not  converge  to  a  steady  growth  rate.  The  growth  rate 
for  FDM-C  is  closest  to  the  consensus.  The  amplitude  profiles  for  scheme  C  are  also 
closest  to  the  consensus  amplitude  profile  for  the  v-component  and  temperature,  and 
schemes  A  and  B  have  much  smaller  magnitudes  than  the  consensus.  The  fmite  differ- 
ence models  are  closer  than  the  finite  element  models  to  the  consensus  for  schemes  A 
and  C  except  for  FEM-C  at  the  top  layer.  FEM-B  is  closer  than  FDM-B.  For  the  v- 
component  phase,  all  of  the  models  depict  the  negative  tilt  of  the  consensus  profile. 

At  2000  km,  the  growth  rates  for  grid  C  are  greater  than  the  consensus,  and 
the  growth  rates  for  grids  A  and  B  are  much  less  than  the  consensus  growth  rate.  The 
growth  rates  for  grid  B  do  not  converge  to  a  steady  growth  rate.  The  growth  rate  for 
FEM-C  is  closest  to  the  consensus.  The  finite  element  models  are  closer  to  the  consensus 
amplitude  profile  for  temperature  and  v-component  than  are  the  fmite  difference  models 
for  grids  B  and  C,  but  not  for  grid  A.  Scheme  C  has  a  greater  magnitude  than  the  con- 
sensus, especially  FDM-C.  FDM-A  is  closest  to  the  consensus  in  magnitude  but  does 
not  have  enough  resolution  to  depict  the  low  level  shear. 
3.     Hyperbolic  Shear  Type  III 

For  this  part  of  the  experiment,  a  =  4.0  and  k  =  1.0.  This  causes  the  strongest 
shear  to  occur  at  about  Z  =  0.25  instead  of  at  the  surface,  and  provides  a  variety  of 
responses  that  are  less  sensitive  to  the  lower  boundary  conditions. 
a.     Sixty-Layer  Models 

At  the  4000  km  wavelength,  all  the  growth  rates  are  identical  except  for 
FEM-B  which  has  a  smaller  growth  rate  (Table  1).  All  of  the  sixty-layer  models  except 
FDM-B  and  FEM-B  converge  to  the  same  solution  for  the  v-component  and  temper- 
ature amplitude  profiles.  FDM-B  is  identical  to  the  consensus  profiles  except  at  the 
surface  of  the  temperature  profile.  FEM-B  shows  the  same  results  as  the  Type  II  case, 
but  the  temperature  oscillations  are  smaller  (Fig.  24).  The  magnitudes  of  the  growth 
rates  and  the  amplitude  profiles  are  larger  than  those  for  Type  I  and  Type  II.  The  phase 
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profiles  of  all  the  models  converge  to  the  same  solution  except  FEM-B  which  has  a 
slightly  larger  phase. 

For  3000  km,  FDM-B,  FDM-C  and  FEM-C  have  identical  growth  rates. 
FDM-A  and  FEM-A  also  have  identical  but  slightly  smaller  growth  rates  and  FEM-B 
has  the  smallest  growth  rate.  The  v-component  and  temperature  amplitude  profiles 
show  small  differences  between  the  models,  especially  near  the  strongest  shear  region  at 
Z  =  0.25  (Figs.  25-26).  For  the  v-component  amplitude  profiles,  the  finite  difierence 
models  converge  to  the  same  solution  and  represent  the  consensus  solution.  FEM-A 
has  a  slightly  smaller  amplitude  and  FEM-C  has  a  slightly  larger  amplitude  in  the  shear 
zone.  In  the  temperature  profiles,  the  differences  are  shghtly  larger  for  FEM-A  and 
FEM-C.  The  results  for  FEM-B  are  the  same  as  the  previous  cases,  but  again  the  os- 
cillations are  smaller.  For  the  v-component  phase  profiles,  all  the  models  converge  to 
the  same  phase  except  FEM-A  and  FEM-B  which  are  slightly  larger  (Fig.  25). 

At  the  2000  km  wavelength,  the  grovMh  rates  show  a  little  more  variation 
than  in  the  other  cases.  FDM-B  and  FDM-C  are  identical  and  represent  the  consensus 
growth  rate.  Although  very^  close  to  the  consensus,  the  growth  rate  for  FDM-A  is 
slightly  smaller  and  FEM-C  is  a  little  larger.  FEM-A  has  a  slightly  smaller  growth  rate 
than  FEM-A  and  FEM-B  has  the  smallest  growth  rate.  The  v-component  and  temper- 
ature amplitudes  also  show  more  variations  among  the  models  in  the  region  of  the 
stongest  shear  than  in  previous  sixty-layer  cases  (Figs  27-28).  FDM-B  and  FDM-C  are 
identical  and  represent  the  consensus  profile.  FDM-A  is  very  close  to  the  consensus, 
but  slightly  smaller  in  magnitude.  FEM-A  has  a  smaller  magnitude  than  the  consensus 
and  FDM-A.  FEM-C  has  a  larger  magnitude  than  the  consensus.  FEM-B  has  a  much 
smaller  magnitude  than  the  consensus.  The  differences  among  the  models  are  slightly 
greater  for  the  temperature  profiles  than  the  v-component  profiles.  However,  for  the 
phase  profiles,  the  finite  difference  models  and  FEM-C  converge  to  the  same  phase. 
FEM-A  is  slightly  greater  and  FEM-B  is  sUghtly  smaller  (Fig.  27). 
b.     Six-Layer  Models 

To  save  space,  the  profiles  for  the  six-  and  two-layer  models  are  not  shown. 
Table  4  contains  the  exponential  growth  rates  the  six-,  four-  and  two-layer  models  for 
Type  III.  For  a  wavelength  of  4000  km,  the  growth  rates  for  the  finite  element  models 
are  closer  to  the  consensus  than  are  the  growth  rates  for  the  fmite  difference  models  for 
schemes  A  and  C,  but  not  B.  FDM-B  is  closest  to  the  consensus  growth  rate.  Overall, 
the  fmite  element  models  are  closer  than  the  fmite  difierence  models  to  the  consensus 
amplitude  profile  for  v-component  and  temperature  for  grids  A  and  C.    For  grid  B,  the 
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fmite  diflerence  model  is  closer  to  the  consensus  than  the  finite  element  model  except 
at  the  lowest  level  above  the  surface.  In  general,  grid  B  is  closest  to  the  consensus  am- 
plitude profile.  For  the  phase  profiles,  all  of  the  models  are  close  to  the  consensus  pro- 
file, but  grid  B  is  closest. 


Table  4 

Exponential  Growth  Rates  (  x  10-^) 

Hyperbolic  Shear  Type  III 

Models 

Consensus 

FDM-A 

FD.M-B 

FDM-C 

FEM-A 

FEM-B 

FEM-C 

6  Layers 

4000  km 
3000  km 
2000  km 

2.89 
3.22 
3.03 

2.56 
2.70 
2.08 

2.80 
2.96 
2.38 

2.60 

2.55 
0.98** 

2.71 
2.99 
2.98 

2.78 
3.02 
2.91 

2.71 

2.89 

1.31** 

4  Layers 

4000  km 
3000  km 
2000  km 

2.89 
3.22 
3.03 

2.26 
1.81 
1.65 

2.93 
3.40 
3.76 

2.48 
2.44 
1.95 

2.48 
2.49 
1.54 

2.93 

3.45 
3.95 

2.48 
2.68 

1.77 

2  Layers 

4000  km 
3000  km 
2000  km 

2.89 

3.22 
3.03 

2.23 
2.55 
2.68 

2.13 
2.13 
1.02-- 

2.53 
3.23 
4.20 

1.61 
1.63 
1.41 

2.37 

2.44 

0.27* 

2.24 
2.76 
3.59 

•"•'  These  cases  do  not  converge  to  a  steady  growth  rate 
**  These  models  are  still  converging  to  a  steady  growth  rate 

For  the  3000  km  case,  the  growth  rates  indicate  that  all  the  models  are  at 
least  8  percent  less  than  the  consensus  growth  rate.  The  fmite  element  models  are  closer 
than  the  fmite  different  models  to  the  consensus  and  FEM-B  is  closest.  Overall,  the  fi- 
nite element  models  are  closer  than  the  finite  difference  models  to  the  consensus  ampli- 
tude profile  for  v-component  and  temperature.  The  only  exceptions  are  at  the  top  layer 
of  the  v-component  amplitude  profiles  and  at  the  midpoint  of  the  temperature  profiles 
where  FDM-B  is  closer.  In  general,  FEM-B  is  closest  to  the  consensus.  For  the  phase 
profiles,  the  A  and  B  profiles  are  closer  than  the  C  profiles  to  the  consensus,  and 
FDM-A  is  closest  to  consensus  phase  profile  for  the  v-component. 
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For  a  wavelength  of  2000  km,  the  growth  rates  for  the  finite  element  models 
are  closer  than  the  growth  rates  for  the  finite  difierence  models  to  the  consensus  growth 
rate  with  FEM-A  being  closest.  The  growth  rates  for  the  finite  difference  models  and 
FEM-C  are  much  less  than  the  growth  rate  for  the  consensus.  The  growth  rates  for  the 
C  models  are  still  converging  at  96  hours,  but  even  after  192  hours  the  growth  rates  in- 
crease only  slightly.  The  amplitude  profiles  for  the  v-component  and  temperature  also 
show  that  the  finite  element  models  are  closer  than  the  finite  difierence  models  to  the 
consensus.  FEM-A  is  closest  to  the  consensus  amplitude  profile.  For  the  phase  profiles, 
the  scheme  A  models  are  closer  to  the  consensus  phase  profile  than  are  the  scheme  B 
and  C  models. 

c.     Four- Layer  Models 

At  an  x-wavelength  of  4000  km,  the  growth  rates  for  FDM-B  and  FEM-B 
are  the  same  and  closest  to  the  consensus  growth  rate.  For  schemes  A  and  C  the  growth 
rates  for  the  finite  element  models  are  closer  to  the  consensus  than  are  the  grovnh  rates 
for  the  finite  difference  models.  The  growth  rates  for  schemes  A  and  C  are  much  less 
than  the  consensus.  The  finite  element  models  are  closer  than  the  finite  difierence 
models  to  the  consensus  amplitude  profile  for  grids  A  and  C  (Fig.  29).  For  grid  B, 
FDM-B  is  closer  than  FEM-B  and  is  the  closest  overall  to  the  consensus  profile.  For 
the  phase  profiles,  all  the  models  follow  the  shape  of  the  consensus  profile  well  (Fig.  29). 
FEM-C  is  closest  to  the  consensus  phase  profile. 

For  3000  km,  the  growth  rates  for  the  finite  element  models  are  closer  than 
the  finite  difierence  models  to  the  consensus  for  schemes  A  and  C,  but  not  B.  FDM-B 
is  closest  to  the  consensus  and  schemes  A  and  C  are  much  less  than  the  consensus.  For 
grids  A  and  C,  the  finite  element  models  are  closer  than  the  finite  difference  models  to 
the  consensus  amplitude  profile  for  the  v-component  and  temperature  (Figs.  30-31).  For 
grid  B,  FDM-B  is  closer  than  FEM-B  to  the  consensus  profile  except  in  the  upper  layers 
of  the  v-component  profile.  Otherwise,  FDM-B  is  closest  to  the  consensus  profile  for 
this  case.   For  the  phase  profiles,  FEM-C  is  closest  to  consensus  phase  profile  (Fig.  30). 

For  2000  km,  the  growth  rates  show  that  none  of  the  four-layer  model 
growth  rates  are  close  to  the  consensus.  The  growth  rates  for  the  finite  difference  models 
are  closer  than  the  finite  element  models  to  the  consensus  growth  rate.  The  growth  rates 
for  the  grid  B  models  are  much  greater  than  the  consensus  and  the  growth  rates  for  the 
grid  A  and  C  models  are  much  less.  The  growth  rate  for  FDM-B  is  closest  to  the  con- 
sensus growth  rate.  The  amplitude  profiles  for  the  v-component  and  temperature  show 
the  same  results  as  the  growth  rates  (Figs.  32-33).    The  phase  profiles  show  that  all  the 
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models  except  FDM-C  have  the  same  shape  as  the  consensus  profile  and  that  scheme 
C  is  closest  to  the  consensus  (Fig.  32). 
d.     Two- Layer  Models 

For  the  4000  km  wavelength,  the  growth  rates  indicate  that  none  of  the 
models  grow  as  rapidly  as  the  consensus  model  for  this  case.  For  grids  A  and  C,  the 
groA^lh  rates  for  the  finite  difierence  models  are  closer  than  the  rates  for  the  finite  ele- 
ment models  to  the  consensus  growth  rate.  FEM-B  is  closer  to  the  consensus  than 
FDM-B,  and  FDM-C  is  closest  to  the  consensus.  The  finite  element  profiles  are  closer 
than  the  finite  difference  profiles  to  the  consensus  for  scheme  B  and  equal  to  or  slightly 
closer  to  the  consensus  for  scheme  C.  FDM-A  is  closer  than  FEM-A.  FEM-B  is  closest 
to  the  consensus  in  this  case  except  at  the  surface  for  the  temperature  profile  where 
FEM-C  is  better.  However,  all  the  surface  predictions  are  considerably  less  than  the 
consensus  prediction.  For  the  phase  profiles,  all  six  models  reflect  the  negative  tilt  of  the 
v-component  even  though  they  have  only  two  levels.  FEM-C  is  closest  to  the  consensus 
phase  profile. 

For  the  3000  km  case,  the  grovi^h  rates  indicate  the  same  results  as  the  4000 
km  case.  The  amplitude  profiles  show  that  except  for  grid  B,  the  finite  difference  models 
are  closer  than  the  finite  element  models  to  the  v-component  and  temperature  amplitude 
consensus  profile.  FDM-C  is  closest  to  the  consensus  profile.  For  the  phase  profiles, 
all  the  models  follow  the  shape  of  the  consensus  phase  profile  except  FDM-C.  FEM-A 
is  closest  to  the  consensus  phase  profile  for  this  case. 

At  the  2000  km  x-wavelength,  the  growth  rates  indicate  that  the  grid  C 
models  have  a  much  greater  growth  rate  than  the  consensus  model  and  that  FEM-A  and 
the  grid  B  models  have  a  much  smaller  growth  rate.  FDM-A  also  has  a  smaller  growth 
rate  than  the  consensus,  but  it  is  closer  to  the  consensus  than  are  the  other  models. 
Overall,  the  finite  element  models  are  closer  than  the  finite  difference  models  to  the 
consensus  amplitude  profile  for  grids  B  and  C,  but  not  grid  A.  FDM-A  is  closest  to  the 
consensus  profile  in  magnitude,  but  does  not  have  enough  resolution  to  depict  the  low 
level  shear.  For  the  v-component  phase  profiles,  the  models  have  a  negative  tilt  similar 
to  the  consensus  profile,  but  do  not  have  enough  resolution  to  depict  the  shape  of  the 
consensus  profile.   FDM-C  is  closest  to  the  consensus  phase  profile. 

C.     USING  FINITE  DIFFERENCE  METHOD  IN  FINITE  ELEMENT  MODEL 

In    this    experiment,    the    finite    difference    equations    and    coefficients    for   the 
thermodynamic  predictive  equations  are  used  in  the  finite  element  model  for  scheme  A 
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instead  of  the  matrix  equations.  The  model  is  run  with  linear  vertical  shear  in  the  ubar 
field  for  sixty-layer  models  and  with  hyperbolic  Type  III  vertical  shear  in  the  ubar  field 
for  six-  and  sixty-layer  models.  Waves  are  generated  in  the  model  using  an  initial  per- 
turbation of  v'  =  5.0  m/s  in  the  cosine  term.  Horizontal  wavelengths  of  4000  and  2000 
km  are  used.  The  results  are  compared  with  the  unmodified  FDM-A  and  FEM-A  pro- 
files. 

1.  Linear  Vertical  Shear  in  the  ubar  Field 

For  both  4000  and  2000  km  wavelengths,  the  temperature  and  v-component 
amplitude  profile  of  the  modified  FEM-A  model  is  identical  to  the  unmodified  FEM-A 
model  (Figs. 34-37).   The  phase  profiles  show  the  same  results. 

2.  Hyperbolic  Type  HI  Shear  in  the  ubar  Field 

a.  Sixty-Layer  Model 

When  the  x-wavelength  is  4000  km,  the  v-component  amplitude  shows  that 
the  modified  sixty-layer  FEM-C  model  has  a  greater  amplitude  than  the  unmodified 
models  throughout  the  entire  atmosphere  (Fig.  38).  The  temperature  amplitude  of  the 
modified  FEM-C  model  is  greater  than  the  amplitudes  of  the  unmodified  models  in  the 
region  of  greatest  vertical  shear  in  the  ubar  field  (Fig.  39).  The  phase  profiles  show  that 
the  phase  for  the  modified  model  is  identical  to  the  unmodified  FEM-A  model  (Fig. 
38-39).  For  the  2000  km  wavelength,  the  results  are  the  same,  but  the  difference  between 
the  amplitudes  for  the  modified  and  unmodified  models  is  greater  (Fig.  40). 

b.  Six-Layer  Model 

For  the  4000  km  wavelength,  the  temperature  and  v-component  amplitudes 
of  the  modified  FEM-A  model  are  nearly  the  same  as  the  profiles  for  FDM-A  (Figs. 
41-42).  The  phase  profiles  show  the  modified  FEM-C  model  is  closer  to  the  unmodified 
FEM-C  model  except  in  the  middle  layers  of  the  temperature  phase  profile  (Figs.  41-42). 

When  the  x-wavelength  is  2000  km,  the  temperature  and  v-component  am- 
plitude profiles  show  that  the  modified  model  lies  between  FDM-A  and  the  unmodified 
FEM-C  but  is  slightly  closer  to  the  finite  difference  model  (Figs.  43-44).  The  phase 
profile  of  the  modified  FEM-C  is  closer  to  the  unmodified  FEM-C  than  to  the  finite 
difference  model. 

D.     BACKWARD  TILTING  INITIAL  VALUES 

For  this  experiment,  the  initial  phase  in  equation  (2.41)  is  set  to  give  a  negative  tilt  by 
letting  S  =  bZ  for  the  vorticity  and  temperature  initial  conditions.  The  experiment  is 
run  for  the  sixty-layer  FDM-C  model  with  linear  ubar  shear  and  with  hyperboHc  Type 
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Ill  ubar  shear.  Each  case  is  run  with  h  =  n  and  h  =  4n  which  are  called  phase  I  and 
phase  II  respectively.  The  x-wavelengths  of  4000,  3000  and  2000  km  are  used  in  each 
case. 

1.     Linear  Vertical  Shear  in  the  ubar  Field 
a.     X-wavelength  of  4000  km 

The  growth  rates  for  this  case  do  not  converge  to  a  steady  growth  rate  al- 
though they  appear  to  be  converging  for  the  3000  and  4000  km  wavelengths  for  phase 
I  (Table  5).  The  phase  profile  for  phase  I  shows  a  negative  tilt  of  about  110  degrees  at 
12  hours  which  becomes  less  negative  by  72  hours  (Fig.  45).  For  phase  II,  the  phase 
profile  has  a  much  stronger  negative  tilt  at  12  hours  which  also  becomes  more  vertical 
with  time. 


Table  5 

Exponential  Growth  Rates  (  x  10--) 

Linear  Shear  (Negative  Initial  Phase) 

Sixty-layer  FDM-C  Model 

X-wavelength  (km) 

4000 

3000 

2000 

Phase  I 

12  hours 

2.51 

2.81 

2.85 

24  hours 

1.50 

1.60 

0.91 

36  hours 

1.48 

1.34 

-0.37 

48  hours 

1.29 

1.00 

-0.23 

60  hours 

1.24 

0.92 

0.76 

72  hours 

1.21 

0.87 

-0.43 

Phase  11 

12  hours 

1.07 

1.35 

1.88 

24  hours 

0.89 

1.82 

3.49 

36  hours 

2.28 

3.20 

2.62 

48  hours 

2.58 

2.37 

0.87 

60  hours 

2.03 

1.64 

-0.59 

72  hours 

1.64 

1.21 

0.01 
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b.  X-wavelength  of  3000  km 

The  growth  rates  show  the  same  general  results  as  the  4000  km  results  ex- 
cept the  growth  rates  are  smaller  for  phase  1  and  slightly  greater  until  36  hours  for  phase 
II.   The  phase  profiles  give  the  same  results  as  the  4000  km  case  (Fig.  46). 

c.  X-wavelength  of  2000km 

In  this  case,  the  growth  rates  show  no  apparent  convergence  to  a  steady 
growth  rate.  The  amplitude  profiles  show  there  is  growth  for  about  24  hours  and  then 
the  profile  oscillates  back  and  forth  without  growing  (Fig.  47).  This  is  expected  because 
the  Eady  solution  is  not  unstable  for  this  wavelength.  Phase  II  has  the  same  pattern 
except  it  increases  in  magnitude  for  48  hours  before  decreasing  shghlty.  For  the  phase 
profiles,  phase  I  has  a  slightly  negative  tilt  at  12  hours  which  becomes  positive  by  72 
hours  (Fig.  48).  Phase  II  has  a  greater  negative  tilt  than  phase  I  at  12  hours,  but  it  is 
less  than  the  tilt  for  the  4000  and  3000  km  cases.  The  tilt  for  phase  II  becomes  positive 
by  48  hours  and  then  becomes  negative  again  at  72  hours. 
2.     Hyperbolic  Type  III  Shear  in  the  ubar  Field 

a.  X-wavelength  of  4000  km 

The  growth  rates  for  the  first  12  hours  are  much  greater  than  the  growth 
rates  for  the  linear  case,  but  converge  quickly  to  a  steady  growth  rate  for  phase  I  (Table 
6).  For  phase  II.  the  early  growth  rates  are  larger  and  converge  slower  than  the  rates 
for  phase  I.  However,  the  growth  rates  for  both  phase  I  and  phase  II  converge  to  the 
same  value.  The  phase  profiles  show  less  negative  tilt  than  the  linear  case  (Fig.  49). 
Phase  II  has  a  larger  negative  tilt  than  phase  I  at  12  hours  in  the  region  of  the  strongest 
vertical  shear  in  the  ubar  field  (Z  =  0.25).  This  negative  tilt  becomes  less  negative  with 
time  but  the  change  is  less  pronounced  than  in  the  linear  case. 

b.  X-wavelength  of  3000  km 

In  this  case,  the  results  are  similar  to  the  results  for  4000  km.  The  growth 
rates  are  greater  than  at  4000  km.  The  phase  results  are  also  the  same  as  the  4000  km 
case  except  for  less  negative  tilt  at  12  hours  (Fig.  50). 
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Table  6 

Exponential  Growth  Rates  (  x  10"-) 

Hyperbolic  Shear  Type  III  (Negative  Initial  Phase^ 

) 

< 

Sixty-layer  FDM-C  Model 

X-wavelength  (km) 

4000 

3000 

2000 

Phase  I 

12  hours 

4.33 

4.50 

3.93 

24  hours 

2.89 

3.25 

3.20 

36  hours 

2.84 

3.23 

3.06 

4S  hours 

2.88 

3.21 

3.03 

60  hours 

2.89 

3.21 

3.04 

72  hours 

2.89 

3.21 

3.03 

Phase  II 

12  hours 

4.44 

5.46 

6.37 

24  hours 

4.24 

4.30 

3.47 

36  hours 

3.12 

3.29 

3.06 

48  hours 

2.90 

3.22 

3.04 

60  hours 

2.88 

3.21 

3.04 

72  hours 

2.89 

3.21 

3.03 

c.     X-wavelength  of  2000  km 

At  2000  km,  the  same  trend  is  still  evident  although  the  gro"wth  rates  de- 
crease slightly  from  the  3000  km  case  after  12  hours.  The  growth  rates  do  not  oscillate 
as  they  did  in  the  linear  case.  The  phase  profiles  show  the  negative  tilt  at  12  hours  is  less 
than  the  previous  cases  (Fig.  51). 

The  results  of  this  experiment  agree  with  the  findings  of  Farrell  (1985)  that 
large  initial  growth  rates  can  occur  in  the  early  time  steps  and  converge  to  the  most 
unstable  growth.  In  the  future,  this  type  of  time  evolution  should  be  used  to  test  the 
various  schemes  for  various  wavelengths. 
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Fig.  4.      60-layer,  4000  km  linear  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.  5.      60-layer,  4000  km  linear  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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LINEAR  SHEAR  (60-L,  3000KM) 
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Fig.  6.      60-layer,  3000  km  linear  shear  experiment  at  96  h  comparing  v-component 
(top)  and  temperature  (bottom)  amplitudes. 
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Fig.  7.      60-layer,  2000  km  linear  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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8.      60-Iaver.  2000  km  linear  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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Fig.  9.      6-layer,  4000  km  linear  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.   10.      6-Iayer,  4000  km  linear  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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UNEAR  SHEAR  (6-L,  3000KM) 
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Fig.  1 1.      6-layer,  3000  km  linear  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 


48 


UNEAE  SHEAR  (6-L,  3000KM) 


Cu  O' 


II  =■ 

U  n 


MODELS 
CONSENSUS 
FDM-A 
FDM-B 

FEM-A 


FEM-B 


'     FEM-C 


0.0  10.0  20.0  30.0  4O.0  SO.O 

TEMPERATURE  AMPLITUDE  (KELVIN) 


Fig.   12.      6-layer.  3000  km  linear  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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UNEAR  SHEAR  (6-L,  2000KM) 
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13.      6-Iayer.  2000  km  linear  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.   14.      6-la\er.  2000  km  linear  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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HYPERBOUC  TANGENT  UBAR  PROFILES 
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Fig.   15.      Hyperbolic  tangent  ubar  profiles. 
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Fig.   16.      -l-layer,  4000  km  Type  I  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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HYPERBOUC  SHEAR  TYPE  I  (4-L,  3000KM) 
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Fig.   17.      4-layer.  3000  km  Type  I  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.   18.      4-layer.  2000  km  Type  I  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.   19.      4-layer.  4000  km  Type  II  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.  20.      4-layer.  4000  km  Type  II  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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H\TERBOLIC  SHEAR  TiTE  II  (4-L,  3000KM) 


MODELS 
CONSENSUS 
FDM-A 
FDM-B 

FEM-A 


FEM-B 


0.0 


2.0  4.0  6.0  8.0  10.0 

r'-COMPON'ENT  AMPUTITDE  (M/S) 


12.0 

•icf 


HYPERBOUC  SHEAR  TYPE  II  (4-L  3000KM) 


MODEI^ 
CONSENSUS 
FDM-A 
.FDM-B^^^^^^ 

FEM-A 


FEM-B 


»     FEM-C 


45        »0       13S      180      225      270      315      3«0 

V-COMPONENT  PHASE  (DEGREES) 


Fig.  21.      4-layer.  3000  km  Type  II  shear  experiment  at  96  h  comparing  v-component 
amplitudes  (top)  and  phases  (bottom). 
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Fig.  22.      4-layer,  3000  km  Type  II  shear  experiment  at  96  h  comparing  temperature 
amplitudes. 
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HYPERBOUC  SHEAR  TYPE  II  (4-L,  2000KM) 
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Fig.  23. 


4-layer,  2000  km  Type  II  shear  experiment  at  96  h  comparing  v-component 
amplitudes. 
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Fig.  24.      60-layer.  -1000  km  Type  III  shear  experiment  at  96  h  comparing  temper- 
ature amplitudes. 
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Fig.  25.       60-layer,  3000  km  Type  III   shear  experiment  at  96  h  comparing  v- 
componeiit  amplitudes  (top)  and  phases  (bottom). 
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mTERBOLIC  SHEAR  TiTE  III  (60-L,  3000KM) 
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Fig.  26. 


60-layer,  3U00  km  Type  III  shear  experiment  at  96  h  comparing:  temper- 
ature amplitudes. 
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HYPERBOLIC  SHEAR  TYPE  UI  (60-L,  2000KM) 
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Fig.  27.       60-layer,  2000  km  Type  III   shear  experiment  at  96  h  comparing  v- 
coniponent  amplitudes  (top)  and  pliases  (bottom). 
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Fig.  28.      60-layer.  2000  km  Type  III  shear  experiment  at  96  h  comparing  temper- 
ature amplitudes. 
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HYPERBOUC  SHEAR  TYPE  lU  (4-L,  4000KM) 
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Fig.  29. 


4-layer,  4000  km  Type   III   shear  e.xperiment   at  96  h  comparing  v- 
component  amplitudes  (top)  and  piiases  (bottom). 
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Fig.  30. 


4-layer,   3000  km  Type   III   shear  experiment  at  96  h  comparing  v- 
coniponent  amplitudes  (top)  and  phases  (bottom). 
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Fig.  31.      4-Iayer,  3000  km  Type  III  shear  experiment  at  96  h  comparing  temper- 
ature amplitudes. 
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Fig.  32.       4-layer,   2000   km   Type   III   shear  experiment  at  96  h  comparing  v- 
component  amplitudes  (top)  and  phases  (bottom). 
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Fig.  33. 


4-layer,  2000  km  Type  III  shear  experiment  at  96  h  comparing  temper- 
ature amplitudes. 
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Fig.  34.      60-layer.  4000  km  linear  thermodynamic  modification  at  96  h  comparing 
v-component  amplitudes  (top)  and  phases  (bottom). 
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Fig.  35.      60-layer,  4000  km  linear  thermodynamic  modification  at  96  h  comparing 
temperature  amplitudes  (top)  and  phases  (bottom). 
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Fig.  36.      60-layer,  2000  km  linear  thermodynamic  modification  at  96  h  comparing 
v-component  amplitudes  (top)  and  phases  (bottom). 


73 


HEIGHT.  Z=  -LN(P/PO) 

0.0    0.1     0.Z     0.3     0.4     0.5     0.6     0.7     0.6     0.9     1.0 

UNEAR  SHEAR  (60-U  2000KM) 

"^]^>                                                     MODELS 
<=s::r''''^                                                 ^     FEM-A 

J^y                                                 0     FEM-A(MOD) 
"^^^ai                                                     "     FDM-A 

^^ 

).o          0.5           1.0           1.5           2.0          as 
TEMPERATURE  AMPLITUDE  (KELVIN) 

LINEAR  SHEAR  (60-U  2000KM) 


-n 


o 
a. 


z 
I 


II 

N! 


^ . ----——-- 

^ 

It!  °^ 

^ 

6- 

6- 

o 

^^ 

y 

{ 

f 

s^^- 

MODELS 
FEM-A 


FEM-A(MOD) 
FDM-A 


45   00   135   180   225   270   315   3«0 

TEMPERATU'RE  PHASE  (DEGREES) 


Fig.  37. 


60-layer,  2000  km  linear  thermodynamic  modification  at  96  h  comparing 
temperature  amplitudes  (top)  and  phases  (bottom). 
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Fig.  38.      60-layer,  4000  kni  Type  III  thermodynamic  modincation  at  96  h  comparing 
v-component  amplitudes  (top)  and  pliases  (bottom). 
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Fig.  39.      60-layer.  4000  km  Type  III  thermodynamic  modification  at  96  h  comparing 
temperature  amplitudes  (top)  and  phases  (bottom). 
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Fig.  40.       60-layer,  2000  km  Type  III  thermodynamic  modification  comparing  v- 
coniponent  (top)  and  temperature  (bottom)  amplitudes. 
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Fig.  41.      6-Iayer,  4000  km  Type  III  thermodynamic  modification  at  96  h  comparing 
v-component  amplitudes  (top)  and  pliases  (bottom). 
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Fig.  42.      6-layer,  4000  km  Type  HI  thermodynamic  modification  at  96  h  comparing 
temperature  amplitudes  (top)  and  phases  (bottom). 
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Fig. 


43.      6-layer,  2000  km  Type  III  thermodynamic  modification  at  96  h  comparing 
v-coniponent  amplitudes  (top)  and  phases  (bottom). 
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Fig.  44.      6-layer,  2000  km  Type  III  thermodynamic  modification  at  96  h  comparing 
temperature  amplitudes  (top)  and  phases  (bottom). 
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Fig.  45. 


60-layer,  4000  km  initial  value  modification  (linear)  comparing  phase  I  (top) 
and  pliase  II  (bottom)  negative  phase  tilt. 
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Fig.  46.      60-Iayer,  3000  km  initial  value  modincation  (linear)  comparing  phase  I  (top) 
and  piiase  II  (bottom)  negative  phase  tilt. 
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Fig.  47.      60-layer,  2000  ktn  initial  value  modification  (linear)  comparing  phase  I  (top) 
and  piiase  II  (bottom)  amplitudes. 
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Fig. 


48.      60-layer.  2000  kni  initial  value  modification  (linear)  comparing  phase  I  (top) 
and  piiase  II  (bottom)  negative  phase  tilt. 
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Fig.  49.      60-layer.  4000  km  initial  value  modification  (Type  III)  comparing  phase  I 
(top)  and  phase  II  (bottom)  negative  phase  tilt. 
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Fig.  50. 


60-layer,  3000  kni  initial  value  modification  (Type  III>  comparing  phase  I 
(top)  and  phase  II  (bottom)  negative  phase  tilt. 
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Fig.  51.      60-layer,  2000  km  initial  value  modification  (Type  III)  comparing  phase  I 
(top)  and  phase  II  (bottom)  negative  phase  tilt. 
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IV.     CONCLUSIONS 

With  a  linear  vertical  ubar  profile,  FEM-B  still  has  an  oscillation  near  the  lower 
boundar>'  of  the  temperature  profile  after  correcting  the  length  of  the  vertical  velocity 
vector  and  is  the  poorest  of  the  six  models.  FEM-A  is  better  than  FEM-B  but  is  still 
not  equal  to  the  consensus  profile.  In  the  hyperbolic  tangent  vertical  ubar  profiles,  the 
finite  element  models  do  not  perform  as  well  as  the  fmite  difference  models  for  Types  I 
and  II  except  for  the  six-layer  model  for  Type  II.  Grid  A  is  closest  to  the  consensus  for 
most  of  the  Type  I  and  Type  II  cases.  However,  for  Type  III,  the  finite  element  models 
are  better  than  the  finite  difference  models  in  many  of  the  cases.  FEM-B  shows  much 
improvement  in  Type  III.  This  may  be  because  the  maximum  vertical  shear  in  the  ubar 
field  occurs  near  Z  =  0.25,  which  may  be  far  enough  from  the  boundar}'  to  enable 
FEM-B  to  overcome  its  earher  problems  near  the  boundar}^ 

Overall,  scheme  C  is  often  better  than  A  and  B,  especially  FEM-C.  Scheme  C  may 
have  an  advantage  because  it  predicts  the  variables  at  one  more  level  than  the  A  and  B 
models  except  for  the  vertical  velocity  and  temperature.  However,  even  for  the  temper- 
ature profile,  where  grid  B  predicts  the  same  number  of  levels  as  grid  C,  scheme  C  often 
has  the  better  surface  prediction.  Grid  A  seems  to  perform  surprisingly  well  in  spite  of 
the  analysis  of  Arakawa  (1984). 

Future  research  should  concentrate  on  the  boundarv'  elements  for  FEM-A  and 
FEM-B.  Perhaps  new  hybrid  elements  could  be  employed.  Another  approach  would 
be  to  continue  the  work  of  isolating  the  effects  of  the  vorticity,  divergence  and 
thermodynamic  terms  in  the  finite  element  models.  This  could  be  accomplished  by  using 
the  finite  difference  method  for  these  equations  in  the  finite  element  models  and  ob- 
serving each  term's  contribution. 

So  far,  all  the  baroclinic  instability  experiments  have  used  a  Tbar  field  that  increases 
hnearly  with  height,  which  may  not  be  too  challenging  for  simple  finite  difference  models 
to  solve.  A  vertical  Tbar  field  that  is  not  linear,  such  as  the  hyperbolic  tangent  profile 
used  for  the  ubar  field  in  this  study,  may  allow  the  finite  element  models  to  take  full 
advantage  of  their  possibiHties.  Finally,  future  work  should  continue  the  initial  value 
growth  experiments  by  evaluating  all  the  models  and  by  using  smaller  resolution  models 
and  smaller  wavelength  in  x.  Smaller  horizontal  wavelengths  may  lead  to  a  different  type 
of  growth. 
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APPENDIX  A.     FINITE  DIFFERENCE  APPROXIMATIONS 

1.    For  terms  of  the  form  —rrr-W: 

dZ 

a.   FDM-A  and  FDM-B,  at  level  Z  =  Z,: 

du  1  ^+1  ~  ^  Uj  -  u,_i 


b.   FDM-C,  at  level  Z  =  Z',: 

HIT  1  ^1+1  ~  ^i 

dZ  2       '^^  Z',^, -Z',  ^     ^  Z',-Z'._,  ^^ 


i+l       -^  1  ^  1       ^  i_i 


2.   For  terms  of  the  form  -:r^W; 

cL 


a.   FDM-A,  at  level  Z  =  Z.: 


rT  1  T  +,  —  T,  Ti  —  T,_i 


b.    FDM-B,  at  level  Z  =  Z,: 


c.   FDM-C,  at  level  Z  =  Z,: 

iIw  =  -Lw[(  ^'^'~^'-^  )  +  (   ^'~^'-'    )] 
az  2      '^^  ZVi-Z'.  ^^  Z'i-Z',-1 
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3.    For  terms  of  the  form  -r^V: 

uZ 


a.    FDM-A,  at  level  Z  =  Z,: 


dTT  1  ^^i+i  ~  ^.  U,  —  U,    , 

dZ  2      '^^  Z,+i-Z,  -"^^  Z,-Z,_i 


b.    FDM-B,  at  level  Z  =  Z',: 


dii  1  ^1-^1  ~  ^1 


dZ  2  '    '+'        '''  Z,+i-Z, 

c.    FDM-C,  at  level  Z  =  Z',: 


dZ         2     '^^  Z'^+i -Z',  ^^^  Z',-Z',_,  ''^ 
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APPENDIX  B.     GALERKIiN  FORM  OF  FEM-A  PROGNOSTIC 

EQUATIONS 

1.   Vorticity  Equations  (2.25)  and  (2.26): 


i+l 


dM  r^ 


i+l 


dt 


j=i-i 


i+l  i+l 


1+1 


^   >    uM   Ay\(/)k<^,dZ-;g  >    V^,j\(/).dZ. 

k=l-I     J=l-1  j=l-l 


(B.l) 


1+1 


1+1 


dAJ 
~di 


(i)4>^CLL  =  — f 


]=i-i 


Zo 


j=l-l 


+ 


i+l  i+l 

■Z-E 

k=i-i     j=i-i 


i+l 


A 


^'4>MA'^-P/  vij\<^idz. 


Zo 


^       *'Zo 

j=i-l 


(B.2) 


Note  that  in  these  equations,  and  the  equations  that  follow,  the  basis  functions  are 
functions  of  Z  ((/),  =  (/),(Z)  and  t/^,  =  (//,(Z))  .  All  of  the  other  variables,  A,  D,  H,  Q,  T,  u, 
U,  V,  and  W,  are  functions  of  time  (A,  =  Ai(t),  D,  =  D,(t),  H,  =  H,(t),  Q,  =  Q,(t), 
Ti  =  T.(t),  u.  =  u,(t),  U.  =  U.(t),  V,  =  V,(t),  W.  =  W,(t))  . 
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2.   Divergence  Equations  (2.27)  and  (2. 28): 


i+l 


VdD'i    f^i 


i+l 


j=i-l 


j=i-i 


(/)(/),dZ 


i+l  1+] 


1+1 


;=i-l      j=i-l  j=i— 1 


i+1  i+l 

k=i-i  j=i-i 


rz 


i+l 


M</>, 


% 


dZ    ^J 


j=i-i 


.A.(/>,dz  +  /  >  h;     (/),(/),dz 


(B.3) 


i+l 


J=l- 


_dDi    rz. 

"dT 


i+l 


'^(^j(/),dZ  =  f\   h\\  '^(/)j</>,dZ 


j=i-i 


i+l  i+l 


i+l 


+  /^ 


y  H*^  y  Di  [^'(/)j0k'/>idZ  -  ^  )    Ui  [^>j<^,dZ 

A^       A_J       ''Zo  Z_j       Jzo 

;=1-1        J=l-1  J=l-1 


i+l  i+l 


rz. 


i+l 


d(/). 


k=l-l         J=l-1  '^O  j=i_l 


M'2.Hij\(/>,( 


(B.4) 
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3.   Thermodynamic  Equations  (2.29)  and  (2.30): 


i+i 


i+l  i+1 


^^JV,'iz=-.^n'^T.j; 


4>i4>k4>A^ 


j=i-i 


k=i-i    ]=i-i 


i+l  i+l 


+ 


k=i-l      j=i-l 


rz. 


% 


d(i). 


dZ    ^J 


(/).(/),dZ 


i+l  i+l 

k=l-l         j=l-l  ^0 


i+l 


dZ 


v!/.(/),dZ  + 


0;(/),dZ. 


(B.5) 


'+1  i+l  i+l 

</)j</),dZ  =  ^ 

]=i-l  k=i-l      j=i-l 


1+1 


1+1 


T^, 


</)  (^k</)idZ 


+ 


i+l  i+l 


R 


V 


d4>, 


dZ    ^J 


</>i0idZ 


k=i-i     ]=i-i 


% 


i+l  i+l 


i+l 


T' 


1^'\ 


^^dc/). 


k=i-i     j=i-i 


Zo 


dZ    ^J 


i/'.^jdZ  + 


EC- 


4>.4>. dZ. 


'j'Ki 


j=i-i 


(B.6) 
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APPENDIX  C.     BASIS  FUNCTION  EQUATIONS  FOR  FEM-A 

1.   Notation: 


A,  =  Z,  -  Z,_i 


A'.  =  Z'i-ZV, 


2.  For  the  general  case: 


(f  -  A',  +  .5A',_,  A',_,  A'i 

-?  +  •5A'.+,  A'j    _  ^  A',^, 


<  c  < 


.5(A',  +  A',^,)  2 

for  any  i  sufficiently  far  from  the  boundaries,  3  <  i  <  n. 


0  <  ^  <  A',+, 


for  2  <  i  <  n. 
3.   For  special  cases,  4)^^  (/j^,  (/>n^i,  </>n-2»  "Ai.  '/'n+i 


-^  +  .5A'2 


A'- 


^+  1.5A'2 


-;  +  •5A'3 
.5(A%  +  A'3) 


-.5AS  <  ^  <  .5A'- 
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'A,(s^)  = 


-c  +  AS 
A', 


0  <  c  <  A'2 


^n+l(0  = 


^  +  A' 


n+l 


A' 


n+l 


<l^n+^iO  = 


.5(A'„^,  +  A',) 


A' 


n+l 


A', 


-A'   <C<- 


A'n+l 


0  <  ^  <  -  A' 


n+l 


0n+2(^)  = 


^  +  .5A' 


n+l 


A' 


n+l 
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APPENDIX  D.     GALERKIN  FORM  OF  FEM-B  PROGNOSTIC 

EQUATIONS 

1.  The  vorticity  equations  have  the  same  form  as  the  vorticity  equations  for  model 
FEM-A,  equations  (B.l)  and  (B.2). 

2.  The  divergence  equations  are  the  same  as  the  divergence  equations  in  model  FEM-A. 
equations  (B.3)  and  (B.4). 

3.  Thermodynamic  Equations  (2.29)  and  (2.30): 


i+l 


ill 

di 


j=i-i 


Zo 


i+l  1+1 

i/'ji/'idZ  =  -  ^  y   u*^  y   T^ 
k=i-i    j=i-i 


•Ajf/'k'AidZ 


Zo 


1+1  i+l 


k=i-l      j=i-l 


'Zt 


'Zo 


dc/), 


dZ   "-^J^' 


4>,^,dZ 


i+l  i+1 


k=i-i    j=i-i 


i+l 


d-At 


% 


dZ    ^J^' 


ij/.ij/AZ  + 


1< 

j=i-i 


ip.ip.dZ. 


(D.l) 
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I+l 


1+1        1+1 


"IT 


^,iPAz  =  ^  ;  u  ;  T 


■jn 


TTk   \    tJ 


!//,(/>  ktA.dZ 


j=i-i 


k=i-i     j=i-i 


+ 


1+1        i+i 


rz 


R 


k=l-l     ]=1-1 


'd<t>y 


•"Zo 


dZ    ^J^' 


(t),ip  ^dZ 


i+l  1+1 

k=i-i     j=i-i 


i+l 


w 


-M«Ai 


dZ  ^J^' 


i/'ji/'.dZ  + 


Z< 


(A.tA.dZ. 


'jri 


J=i- 


(D.2) 


Note  that  in  these  equations  the  basis  functions  are  functions  of  Z  (</>,  =  (/),(Z)  and 
ijy,  =  i/',(Z))  .  All  of  the  other  variables,  Q,  T,  u,  U,  V,  and  W,  are  functions  of  time 
{Q,  =  Q,(t),  T,  =  T.(t),  u,  =  u,(t),  U,  =  U,(t),  V,  =  V,(t),  W,  =  W.(t))  . 
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APPENDIX  E.     GALERKIN  FORM  OF  FEM-C  PROGNOSTIC 

EQUATIONS 

1.  The  vorticity  equations  have  the  same  form  as  the  vorticity  equations  for  model 
FEM-A,  equations  (B.l)  and  (B.2). 

2.  Divergence  Equations  (2.27)  and  (2.28): 

i+l  1+1 


7   -^j\<^.dZ  =  fN   A^,j%j(/),dZ 

J=l-1  j=l-l 


i+I  1+1  i+1 


k=i-i     ]=i-i  j=i-i 


i+l  i+1  i+1 


M 


y  u'^  y  \w{  '-^  (/>j(/),dz + A^'  y  h^ j ^^^^.dz.         (e.  i ) 


7  ^—  ^ 

k=i-i     j=i-i        ^  i=i-i 
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i+1 


ZdD'.    C^ 


i+1 


0,(/),dZ  =  f  >   A\ 


(/).(/)jdZ 


j=i-i 


j=i-i 


i+1  i+1 


i+1 


+  M 


y  u*^  y  Di  j^'c/)j0k0.dz  -  ^  y  ui  j%j0.dz 

;=i— 1      J=i— 1  j=i— 1 


+  )" 


i+1  i+1 

I'-l 

k=l-l     J=l-1 


w 


i+1 


0,0,dZ  +  M'  >    HiJ     c/>j</>,dZ. 


dZ    ^J^' 


J=l-1 


(E.2) 


Note  that  in  these  equations,  and  the  equations  that  follow,  the  basis  functions  are 
functions  of  Z  (0,  =  4>,{Z)  and  )//,  =  ^,{Z))  .  All  of  the  other  variables.  A,  D,  H,  Q,  T,  u, 
U,  V  and  \V,  are  functions  of  time  (A,  =  A,(t),  D,  =  Di(t),  H,  =  H,(t),  Q,  =  Q,(t), 
T,  =  T,(t),  u,  =  u,(t),  U,  =  U.(t),  V,  =  V,(t)  and  W,  =  W,(t))  . 

3.   Thermodynamic  Equations  (2.29)  and  (2.30): 


i+1 


j=i-i 


i+1  1+1 


dt 


'l''l< 


(/),(/)  i((/),dZ 


k=i-i    j=i-i 


i+I  1+1 


+ 


k=i-l      ]=i-l 


r^T 


Zo 


d4>y 


dZ    ^J^' 


(/).</),dZ 


i+1  i+1 


S"Z 


Zt 


i+1 


Wi       -^  <^j</.,dZ  +  y  Ql  j%j0,dZ. 
k=i-i     j=i-i      ^^  i=i-i 


(E.3) 
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1+1 


"dT 


j=i-i 


rz. 


^ 


(/),<f),dZ  =  M 


'JY    1 


i+1  1+1 

k=i-i     j=i-i 


rz. 


T^, 


0,(/)^.0,dZ 


i+l  i+1 


66, 


dZ   ^J^' 


(f).4>, dZ 


i+l  i+l 


i+l 


T' 


k=i-i     ]=i-i 


Jw. 


^"d(^. 


Zo 


dZ    ^J 


0,0idZ  + 


2^Qij\</>idZ. 

j=i-i 


(E.4) 
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